#############################################################
##---------------------------------------------------------##
#############################################################
####THE ORIGINS OF ECONOMIC VOTING: Investigatiing Where and When Vote Swings Come From
####Soledad Artiz and Jonathan Phillips
####5-5-12
####soledadartiz@fas.harvard.edu; jonathanphillips@fas.harvard.edu
####
####Data Provided by Tandon (2012), "Economic reform, voting, and local political intervention: Evidence from India"
####Access to the data was kindly provided by Dr. Sharad Tandon. While we were unable to obtain permission to make the data publicly available, those wishing to replicate our results are kindly requested to contact Dr. Tandon on standon@ers.usda.gov to request access to the data.
#############################################################
##---------------------------------------------------------##
#############################################################

data<-read.csv("BaselineFinal.csv")
names(data)

##Re-organize the data by purpose/type
VSchange<-data[,1:11]
Reforms<-data[,c(1:5, 12:19)]
Falsif<-data[,c(1:5, 19:21)]
Regcontrols<-data[1:118, c(1:5, 22:27)]
 
comp <- cbind(VSchange,Reforms[,-(1:5)])
head(comp)
attach(comp)

vote.changes<-cbind(comp$inc, comp$bjp, comp$com, comp$regional)
colnames(vote.changes)<-c("INC", "BJP","COM","REG")
summary(vote.changes)


#############################################################
##--------- Functions for Clustered Standard Errors--------##
#############################################################
clx <- function(fm, dfcw, cluster){
  library(sandwich)
  library(lmtest)
  M <- length(unique(cluster))
  N <- length(cluster)
  dfc <- (M/(M-1))*((N-1)/(N-fm$rank))
  u <- apply(estfun(fm),2,
             function(x) tapply(x, cluster, sum))
  vcovCL <- dfc*sandwich(fm, meat=crossprod(u)/N)*dfcw
  coeftest(fm, vcovCL) 
}


#############################################################
##----------------------Section 2--------------------------##
##--------------Replicating Tandon Table 4-----------------##
#############################################################

##TABLE 1: Tandon (2012) TABLE 4

##Load PLM package designed to work with panel data
library(plm)

##Specify Length of the Cluster Variable, i.e. how many clusters there are
m<-length(unique(data$region))

##Create Dataframe Indexing by the Panel Variables - region and time
data.plm<-pdata.frame(data, index=c("region", "time"))
##Region is group variable and time is time variable

##Model 1
lm4.1<-plm(incum_voteshare~tarmeasure + factor(time), model="pooling", data=data.plm)
lm4.1
se4.1<-vcovHC(lm4.1, method="arellano", type="HC1", cluster="group")*(m/(m-1))
se4.1<-sqrt(diag(se4.1))
se4.1

##Model 2
lm4.2<-plm(incum_voteshare~tarmeasure + factor(time) + fdi_reform + ps_reform + qr_reform + licensing_reform + tardecreases, model="within", data=data.plm)
lm4.2
se4.2<-vcovHC(lm4.2, method="arellano", type="HC1", cluster="group")*(m/(m-1))
se4.2<-sqrt(diag(se4.2))
se4.2

##Model 3
library(AER)
iv4.3<-ivreg(incum_voteshare~tarmeasure + fdi_reform + ps_reform + qr_reform + licensing_reform + tardecreases + factor(region) + factor(time) | tarmeasure_instrument + factor(time) + fdi_reform + ps_reform + qr_reform + licensing_reform + tardecreases + factor(region) + factor(time), model=TRUE, data=data.plm)
summary(iv4.3)

iv4.3b <- lm4.2 #fake model to allow IV to enter inside apsrtable
iv4.3b$coefficients[1] <- iv4.3$coefficients[2]
iv4.3b$vcov[1,1] <- vcov(iv4.3)[2,2]


##Model 4
lm4.4<-plm(inc~tarmeasure + fdi_reform + licensing_reform + ps_reform + qr_reform + tardecreases + factor(time), model="within", data=data.plm)
lm4.4
vcov4.4<-vcovHC(lm4.4, method="arellano", type="HC1", cluster="group")*(m/(m-1))
se4.4<-sqrt(diag(vcov4.4))
se4.4

llm4.4<-lm(inc~tarmeasure + factor(time) + fdi_reform + ps_reform + qr_reform + licensing_reform + tardecreases + factor(region) , data=data.plm)
llm4.4$se <- c(se4.4,rep(0,length(unique(data.plm$region))))

##Model 5
lm4.5<-plm(com~tarmeasure + fdi_reform + licensing_reform + ps_reform + qr_reform + tardecreases + factor(time), model="within", data=data.plm)
lm4.5
vcov4.5<-vcovHC(lm4.5, method="arellano", type="HC1", cluster="group")*(m/(m-1))
se4.5<-sqrt(diag(vcov4.5))
se4.5

llm4.5<-lm(com~tarmeasure + factor(time) + fdi_reform + ps_reform + qr_reform + licensing_reform + tardecreases + factor(region) , data=data.plm)
llm4.5$se <- c(se4.5,rep(0,length(unique(data.plm$region))))


##Model 6
lm4.6<-plm(bjp~tarmeasure + fdi_reform + licensing_reform + ps_reform + qr_reform + tardecreases + factor(time), model="within", data=data.plm)
lm4.6
vcov4.6<-vcovHC(lm4.6, method="arellano", type="HC1", cluster="group")*(m/(m-1))
se4.6<-sqrt(diag(vcov4.6))
se4.6

llm4.6<-lm(bjp~tarmeasure + factor(time) + fdi_reform + ps_reform + qr_reform + licensing_reform + tardecreases + factor(region) , data=data.plm)
llm4.6$se <- c(se4.6,rep(0,length(unique(data.plm$region))))


##Model 7
lm4.7<-plm(regional~tarmeasure + fdi_reform + licensing_reform + ps_reform + qr_reform + tardecreases + factor(time), model="within", data=data.plm)
lm4.7
vcov4.7<-vcovHC(lm4.7, method="arellano", type="HC1", cluster="group")*(m/(m-1))
se4.7<-sqrt(diag(vcov4.7))
se4.7

llm4.7<-lm(regional~tarmeasure + factor(time) + fdi_reform + ps_reform + qr_reform + licensing_reform + tardecreases + factor(region) , data=data.plm)
llm4.7$se <- c(se4.7,rep(0,length(unique(data.plm$region))))


#############################################################

##Calculating Expected Values

##Set All Covariates at their median or mean
median_tandon<-apply(cbind(data.plm[,13], data.plm[,16:19], data.plm[,15],data.plm[,4:6]),2,FUN=median)

##INC
##Simulate Betas for INC Model
library(mvtnorm)
set.seed(1234)
simbetas_inc <- rmvnorm(1000,mean=lm4.4$coefficients,sigma=vcov4.4)
##Estimate Expected Values
ev.inc<-simbetas_inc %*% median_tandon
summary(ev.inc)

##COM
set.seed(1234)
simbetas_com <- rmvnorm(1000,mean=lm4.5$coefficients,sigma=vcov4.5)
##Estimate Expected Values
ev.com<-simbetas_com %*% median_tandon
summary(ev.com)

##BJP
set.seed(1234)
simbetas_bjp <- rmvnorm(1000,mean=lm4.6$coefficients,sigma=vcov4.6)
##Estimate Expected Values
ev.bjp<-simbetas_bjp %*% median_tandon
summary(ev.bjp)

##Regional
set.seed(1234)
simbetas_reg <- rmvnorm(1000,mean=lm4.7$coefficients,sigma=vcov4.7)
##Estimate Expected Values
ev.reg<-simbetas_reg %*% median_tandon
summary(ev.reg)

##Matrix of Expected Values
expected.values_tandon<-cbind(ev.inc, ev.bjp, ev.com, ev.reg)
colnames(expected.values_tandon)<-c("EV[INC]", "EV[BJP]", "EV[COM]", "EV[REG]")
summary(expected.values_tandon)


##FIGURE NOT INCLUDED: Plot Expected Values
plot(density(expected.values_tandon[,3]), col="gold2", xlim=c(-.02, .02), lwd=2, xlab="Simulated Expected Values", sub="Covariates at Medians", main="Expected Values using Table 4 \n Tandon (2012)")
lines(density(expected.values_tandon[,1]), col="magenta4", lwd=2)
lines(density(expected.values_tandon[,2]), col="darkseagreen4", lwd=2)
lines(density(expected.values_tandon[,4]), col="orangered3", lwd=2)
par(font=2)
legend(x="topright", legend=c("BJP","INC", "Far Left", "Regional"), col=c("darkseagreen4", "magenta4", "gold2", "orangered3"), lty=c(1,1,1,1), cex=c(.7,.7,.7,.7), lwd=c(2,2,2,2))

#############################################################

##Calculating First Differences

#Set values of Tariff Measue at 0 and minimum
notar_tandon <- median_tandon
notar_tandon[2] <- 0
tarmin_tandon <- median_tandon
tarmin_tandon[2] <- min(tarmeasure)

#INC
fd.inc<-(simbetas_inc %*% notar_tandon) - (simbetas_inc %*% tarmin_tandon)
summary(fd.inc)

#COM
fd.com<-(simbetas_com %*% notar_tandon) - (simbetas_com %*% tarmin_tandon)
summary(fd.com)

##BJP
fd.bjp<-(simbetas_bjp %*% notar_tandon) - (simbetas_bjp %*% tarmin_tandon)
summary(fd.bjp)

#REG
fd.reg<-(simbetas_reg %*% notar_tandon) - (simbetas_reg %*% tarmin_tandon)
summary(fd.reg)

##Create Vector of First Differences
first.differences_tandon<-cbind(fd.inc, fd.bjp, fd.com, fd.reg)
colnames(first.differences_tandon)<-c("FD[INC]", "FD[BJP]", "FD[COM]", "FD[REG]")
summary(first.differences_tandon)

##FIGURE NOT INCLUDED: Plot First Differences
plot(density(first.differences_tandon[,3]), col="gold2", xlim=c(-0.2, 0.2), lwd=2, xlab="Simulated First Difference", sub="No Tariff - Min Tariff", main="First Differences using Table 4 \n Tandon (2012)")
lines(density(first.differences_tandon[,1]), col="magenta4", lwd=2)
lines(density(first.differences_tandon[,2]), col="darkseagreen4", lwd=2)
lines(density(first.differences_tandon[,4]), col="orangered3", lwd=2)
par(font=2)
legend(x="topright", legend=c("BJP","INC", "Far Left", "Regional"), col=c("darkseagreen4", "magenta4", "gold2", "orangered3"), lty=c(1,1,1,1), cex=c(.7,.7,.7,.7), lwd=c(2,2,2,2))


#############################################################
##----------------------Section 3--------------------------##
##--------Accounting for Dynamics: Including LDV ----------##
#############################################################

##TABLE 2
##Including a Lagged Dependent Variable into the key models from table 4 (see replication file for code for these models)

##Generate Lags
data.plm$incum_1<-lag(data.plm$incum_voteshare,1)
data.plm$left_1<-lag(data.plm$left,1)
data.plm$bjp_1<-lag(data.plm$bjp,1)
data.plm$inc_1<-lag(data.plm$inc,1)
data.plm$com_1<-lag(data.plm$com,1)
data.plm$regional_1<-lag(data.plm$regional,1)


##Model 1
##Using his exact specification and adding LDV: generated LDV (excludes one observation)
lm4.1.ldv<-plm(incum_voteshare~tarmeasure + incum_1 + factor(time), model="pooling", data=data.plm)
se4.1.ldv<-sqrt(diag(vcovHC(lm4.1.ldv, method="arellano", type="HC1", cluster="group")*(m/(m-1))))
summary(lm4.1.ldv)
se4.1.ldv

##Model 2
##Using his exact specification and adding LDV: generated LDV (excludes one observation)
lm4.2.ldv<-plm(incum_voteshare~tarmeasure + incum_1 + factor(time) + fdi_reform + ps_reform + qr_reform + licensing_reform + tardecreases, model="within", data=data.plm)
se4.2.ldv<-sqrt(diag(vcovHC(lm4.2.ldv, method="arellano", type="HC1", cluster="group")*(m/(m-1))))
summary(lm4.2.ldv)
se4.2.ldv

##Model 4
##Using his exact specification and adding LDV: generated LDV (excludes one observation)
lm4.4.ldv<-plm(inc~tarmeasure + inc_1 + factor(time) + fdi_reform + ps_reform + qr_reform + licensing_reform + tardecreases, model="within", data=data.plm)
se4.4.ldv<-sqrt(diag(vcovHC(lm4.4.ldv, method="arellano", type="HC1", cluster="group")*(m/(m-1))))
summary(lm4.4.ldv)
se4.4.ldv

##Model 5
##Using his exact specification and adding LDV: generated LDV (excludes one observation)
lm4.5.ldv<-plm(com~tarmeasure + com_1 + factor(time) + fdi_reform + ps_reform + qr_reform + licensing_reform + tardecreases, model="within", data=data.plm)
se4.5.ldv<-sqrt(diag(vcovHC(lm4.5.ldv, method="arellano", type="HC1", cluster="group")*(m/(m-1))))
summary(lm4.5.ldv) ##SIGNIFICANCE DISAPPEARS!
se4.5.ldv

##Model 6
##Using his exact specification and adding LDV: generated LDV (excludes one observation)
lm4.6.ldv<-plm(bjp~tarmeasure + bjp_1 + factor(time) + fdi_reform + ps_reform + qr_reform + licensing_reform + tardecreases, model="within", data=data.plm)
se4.6.ldv<-sqrt(diag(vcovHC(lm4.6.ldv, method="arellano", type="HC1", cluster="group")*(m/(m-1))))
summary(lm4.6.ldv) ##SIGNIFICANCE DISAPPEARS!
se4.6.ldv

##Model 7
##Using his exact specification and adding LDV: generated LDV (excludes one observation)
lm4.7.ldv<-plm(regional~tarmeasure + regional_1 + factor(time) + fdi_reform + ps_reform + qr_reform + licensing_reform + tardecreases, model="within", data=data.plm)
se4.7.ldv<-sqrt(diag(vcovHC(lm4.7.ldv, method="arellano", type="HC1", cluster="group")*(m/(m-1))))
summary(lm4.7.ldv) 
se4.7.ldv

#############################################################
##----------------------Section 4--------------------------##
##---------------Composition Voting Models-----------------##
#############################################################

##FIGURE 1: Table Showing Overprediction/Underprediction by Tandon Model
VS <- llm4.4$fitted.values + llm4.5$fitted.values + llm4.6$fitted.values + llm4.7$fitted.values
summary(VS)
par(mfrow=c(1,1))
pdf(file="Overpredict.pdf", width = 5, height = 5, family = "Helvetica", pointsize = 10)
plot(VS,pch=3,cex=0.75,col="green",xlab="Observations",ylab="Sum of Change in Predicted Vote Shares") 
abline(a=0,b=0)
dev.off()

## Transforming the data to range from 0 to 1
inc_tran4<-(comp$inc + 1)/4
bjp_tran4<-(comp$bjp + 1)/4
com_tran4<-(comp$com + 1)/4
reg_tran4<-(comp$regional + 1)/4
test_sum4<-inc_tran4 + bjp_tran4 + com_tran4 + reg_tran4
summary(test_sum4)

##Transforming the data using a multivariate logit
bjp_likeli2<-log(bjp_tran4/inc_tran4)
com_likeli2<-log(com_tran4/inc_tran4)
reg_likeli2<-log(reg_tran4/inc_tran4)
summary(bjp_likeli2)
summary(com_likeli2)
summary(reg_likeli2)

##TEST UNTRANSFORMING TRANSFORMATION 
vbjp_tran4<-exp(bjp_likeli2)/(1 + exp(bjp_likeli2) + exp(com_likeli2) + exp(reg_likeli2))
summary(vbjp_tran4)  ##Gets us back to bjp_tran4!
vote_bjp_tran4<-vbjp_tran4*4-1
summary(vote_bjp_tran4) ## Slightly off but gets us back to comp$bjp

##Test 2
vote_bjp2_tran4<-(4*comp$bjp + 4)/(comp$inc + comp$bjp + comp$com + comp$regional + 4)-1
summary(vote_bjp2_tran4)  ##GOOD - same as vote_bjp!!

##Test 3
vote_bjp3_tran4<-((bjp_tran4)/(inc_tran4 + bjp_tran4 + com_tran4 + reg_tran4))*4 -1
summary(vote_bjp3_tran4)  ##GOOD - same as vote_bjp!!


#############################################################
##---------Compositional Voting Analysis-------------------##
#############################################################

##Actual Analysis

#Run separate regressions using OLS
m1_bjp <- lm(bjp_likeli2 ~ tarmeasure + fdi_reform + licensing_reform + ps_reform + qr_reform + tardecreases + factor(time) + factor(region))
m1_com <- lm(com_likeli2 ~ tarmeasure + fdi_reform + licensing_reform + ps_reform + qr_reform + tardecreases + factor(time) + factor(region))
m1_reg <- lm(reg_likeli2 ~ tarmeasure + fdi_reform + licensing_reform + ps_reform + qr_reform + tardecreases + factor(time) + factor(region))
summary(m1_bjp)
#Coefficient in m1R describes how log-ratio of BJP vote share changes relative to INC. Unintuitive.
#Beta QOI is a complicated function of original log transformation (King and Katz). Hard analytically, therefore simulate.


#############################################################

#SUR model
library(systemfit)
one <- bjp_likeli2 ~ tarmeasure + fdi_reform + licensing_reform + ps_reform + qr_reform + tardecreases + factor(time) + factor(region)
two <- com_likeli2 ~ tarmeasure + fdi_reform + licensing_reform + ps_reform + qr_reform + tardecreases + factor(time) + factor(region)
three <- reg_likeli2 ~ tarmeasure + fdi_reform + licensing_reform + ps_reform + qr_reform + tardecreases + factor(time) + factor(region)
form <- list(one, two, three)
SURM <- systemfit(form,"SUR")
summary(SURM)
names(SURM)
vcov.SURM<-vcov(SURM)

#compare SUR and OLS models
rbind(m1_bjp$coefficients,m1_com$coefficients,m1_reg$coefficients) #absolutely identical
vcov(m1_bjp)
SURM$coefCov #Covariance matrix within each regression also seems exactly the same. Of course we do now have additional entries, since vcov is 9 times as big...


#############################################################
##-----------Simulating Quantities of Interest-------------##
#############################################################

##1. Expected Values

##Set All Covariates at their median or mean
median<-c(1,apply(cbind(comp[,13],comp[,16:19], comp[,15], comp[,4:6]),2,FUN=median), rep(0,117))

##BJP
##Simulate Betas for BJP Model
library(mvtnorm)
set.seed(1234)
simbetasm1_bjp <- rmvnorm(1000,mean=SURM$coefficients[1:127],sigma=vcov.SURM[1:127,1:127])

##Estimate Expected Values
ev.m1_bjp<-simbetasm1_bjp %*% median
summary(ev.m1_bjp)


##COM
##Simulate Betas for COM Model
set.seed(1234)
simbetasm1_com <- rmvnorm(1000,mean=SURM$coefficients[128:254],sigma=vcov.SURM[128:254,128:254])

##Estimate Expected Values
ev.m1_com<-simbetasm1_com %*% median
summary(ev.m1_com)


##REGIONAL
##Simulate Betas for Regional Model
set.seed(1234)
simbetasm1_reg <- rmvnorm(1000,mean=SURM$coefficients[255:381],sigma=vcov.SURM[255:381,255:381])

##Estimate Expected Values
ev.m1_reg<-simbetasm1_reg %*% median
summary(ev.m1_reg)


#Untransform
untran_bjp <- exp(ev.m1_bjp)/(1+exp(ev.m1_bjp)+exp(ev.m1_com)+exp(ev.m1_reg))
untran_com <- exp(ev.m1_com)/(1+exp(ev.m1_bjp)+exp(ev.m1_com)+exp(ev.m1_reg))
untran_reg <- exp(ev.m1_reg)/(1+exp(ev.m1_bjp)+exp(ev.m1_com)+exp(ev.m1_reg))
untran_inc <- 1-(untran_bjp + untran_com + untran_reg)

exp_bjp <- (untran_bjp*4)-1
exp_com <- (untran_com*4)-1
exp_reg <- (untran_reg*4)-1
exp_inc <- (untran_inc*4)-1

expected.values<-cbind(exp_inc, exp_bjp, exp_com, exp_reg)
colnames(expected.values)<-c("EV[INC]", "EV[BJP]", "EV[COM]", "EV[REG]")
summary(expected.values)
summary(vote.changes)


##FIGURE NOT INCLUDED: Plot Expected Values
pdf("ev1.pdf")
plot(density(expected.values[,3]), col="gold2", xlim=c(-0.5, 0.5), lwd=2, xlab="Simulated Expected Values", sub="Covariates at Medians", main="Expected Values using Compositional \n Voting Models")
lines(density(expected.values[,1]), col="magenta4", lwd=2)
lines(density(expected.values[,2]), col="darkseagreen4", lwd=2)
lines(density(expected.values[,4]), col="orangered3", lwd=2)
par(font=2)
legend(x="topright", legend=c("BJP","INC", "Far Left", "Regional"), col=c("darkseagreen4", "magenta4", "gold2", "orangered3"), lty=c(1,1,1,1), lwd=c(2,2,2,2))
dev.off()


##FIGURE 2: Plot Side by Side with Tandon's Expected Values
pdf("evandtandonev.pdf")
par(mfrow=c(1,2))
##Expected Values from Tandon's Table 4
plot(density(expected.values_tandon[,3]), col="gold2", xlim=c(-.05, .05), lwd=2, xlab="Simulated Expected Values", sub="Covariates at Medians", main="Expected Values using Table 4 \n Tandon (2012)")
lines(density(expected.values_tandon[,1]), col="magenta4", lwd=2)
lines(density(expected.values_tandon[,2]), col="darkseagreen4", lwd=2)
lines(density(expected.values_tandon[,4]), col="orangered3", lwd=2)
##Expected Values from SUR
plot(density(expected.values[,3]), col="gold2", xlim=c(-0.5, 0.5), lwd=2, xlab="Simulated Expected Values", sub="Covariates at Medians", main="Expected Values using Compositional \n Voting Models")
lines(density(expected.values[,1]), col="magenta4", lwd=2)
lines(density(expected.values[,2]), col="darkseagreen4", lwd=2)
lines(density(expected.values[,4]), col="orangered3", lwd=2)
par(font=2)
legend(x="topright", legend=c("BJP","INC", "Far Left", "Regional"), col=c("darkseagreen4", "magenta4", "gold2", "orangered3"), lty=c(1,1,1,1), lwd=c(2,2,2,2), cex=c(.45,.45,.45,.45))
dev.off()


#############################################################

##2. FIRST DIFFERENCES

#Set values of Tariff Measue at 0 and minimum
notar <- median
notar[2] <- 0
tarmin <- median
tarmin[2] <- min(tarmeasure)

#Calculate Expected Values Given Certain Levels of Tariff Reform
notarcom <- (simbetasm1_com %*% notar)
tarcom <- (simbetasm1_com %*% tarmin)
notarreg <- (simbetasm1_reg %*% notar)
tarreg <- (simbetasm1_reg %*% tarmin)
notarbjp <- (simbetasm1_bjp %*% notar)
tarbjp <- (simbetasm1_bjp %*% tarmin)

##Untransform Expected Values
untran_notarcom <- exp(notarcom)/(1+exp(notarbjp)+exp(notarcom)+exp(notarreg))
untran_notarreg <- exp(notarreg)/(1+exp(notarbjp)+exp(notarcom)+exp(notarreg))
untran_notarbjp <- exp(notarbjp)/(1+exp(notarbjp)+exp(notarcom)+exp(notarreg))
untran_notarinc <- 1-(untran_notarbjp + untran_notarcom + untran_notarreg)

notar_fd_bjp <- (untran_notarbjp*4)-1
notar_fd_com <- (untran_notarcom*4)-1
notar_fd_reg <- (untran_notarreg*4)-1
notar_fd_inc <- (untran_notarinc*4)-1

untran_tarcom <- exp(tarcom)/(1+exp(tarbjp)+exp(tarcom)+exp(tarreg))
untran_tarreg <- exp(tarreg)/(1+exp(tarbjp)+exp(tarcom)+exp(tarreg))
untran_tarbjp <- exp(tarbjp)/(1+exp(tarbjp)+exp(tarcom)+exp(tarreg))
untran_tarinc <- 1-(untran_tarbjp + untran_tarcom + untran_tarreg)

tar_fd_bjp <- (untran_tarbjp*4)-1
tar_fd_com <- (untran_tarcom*4)-1
tar_fd_reg <- (untran_tarreg*4)-1
tar_fd_inc <- (untran_tarinc*4)-1

##Calculate First Differences from Untransformed Values
fd2_bjp <- notar_fd_bjp - tar_fd_bjp
fd2_com <- notar_fd_com - tar_fd_com
fd2_reg <- notar_fd_reg - tar_fd_reg
fd2_inc <- notar_fd_inc - tar_fd_inc

first.differences<-cbind(fd2_inc, fd2_bjp, fd2_com, fd2_reg)

colnames(first.differences)<-c("FD[INC]", "FD[BJP]", "FD[COM]", "FD[REG]")
summary(first.differences)
summary(vote.changes)


##FIGURE NOT INCLUDED: Plot First Differences
pdf("fd1.pdf")
par(mfrow=c(1,1))
plot(density(first.differences[,3]), col="gold2", xlim=c(-0.4, 0.4), lwd=2, xlab="Simulated First Difference", sub="EV[Vote Change | No Tariff] - EV[Vote Change | Min Tariff]", main="First Differences using Compositional \n Voting Models")
lines(density(first.differences[,1]), col="magenta4", lwd=2)
lines(density(first.differences[,2]), col="darkseagreen4", lwd=2)
lines(density(first.differences[,4]), col="orangered3", lwd=2)
par(font=2)
legend(x="topright", legend=c("BJP","INC", "Far Left", "Regional"), col=c("darkseagreen4", "magenta4", "gold2", "orangered3"), lty=c(1,1,1,1), lwd=c(2,2,2,2))
dev.off()

##FIGURE NOT INCLUDED: Plot Side by Side with Tandon's First Differences
pdf("fdandtandonfd.pdf")
par(mfrow=c(1,2))
##Tandon's First Differences
plot(density(first.differences_tandon[,3]), col="gold2", xlim=c(-0.4, 0.4), lwd=2, xlab="Simulated First Difference", sub="No Tariff - Min Tariff", main="First Differences using Table 4 \n Tandon (2012)")
lines(density(first.differences_tandon[,1]), col="magenta4", lwd=2)
lines(density(first.differences_tandon[,2]), col="darkseagreen4", lwd=2)
lines(density(first.differences_tandon[,4]), col="orangered3", lwd=2)
##SUR First Differences
plot(density(first.differences[,3]), col="gold2", xlim=c(-0.4, 0.4), lwd=2, xlab="Simulated First Difference", sub="No Tariff - Min Tariff", main="First Differences using Compositional \n Voting Models")
lines(density(first.differences[,1]), col="magenta4", lwd=2)
lines(density(first.differences[,2]), col="darkseagreen4", lwd=2)
lines(density(first.differences[,4]), col="orangered3", lwd=2)
par(font=2)
legend(x="topright", legend=c("BJP","INC", "Far Left", "Regional"), col=c("darkseagreen4", "magenta4", "gold2", "orangered3"), lty=c(1,1,1,1), cex=c(.5,.5,.5,.5), lwd=c(2,2,2,2))
dev.off()


#############################################################
##--------------Plotting Hexagonal Simplexes---------------##
#############################################################

#A function for transforming between (a,b,c) and (x,y) for vectors
terncoord <- function(a){
  y <- a[,1]*(sqrt(3)/2)
  x <- y*tan(30*pi/180)+a[,3]*sqrt(3)/(2*sin(60*pi/180))
  out <- cbind(x,y)
  return(out)
}

#Function for scalars
terncoords <- function(a){
  y <- a[1]*(sqrt(3)/2)
  x <- y*tan(30*pi/180)+a[3]*sqrt(3)/(2*sin(60*pi/180))
  out <- cbind(x,y)
  return(out)
}


###############################################################

##Figure 3: Basic Hexagon of Predicted Values

#Predicted values in original scale
fv_bjp <- m1_bjp$fitted.values
fv_com <- m1_com$fitted.values
fv_reg <- m1_reg$fitted.values
fv_inc <- 0 - (m1_bjp$fitted.values + m1_com$fitted.values + m1_reg$fitted.values)

untran_fv_bjp <- exp(fv_bjp)/(1+exp(fv_bjp)+exp(fv_com)+exp(fv_reg))
untran_fv_com <- exp(fv_com)/(1+exp(fv_bjp)+exp(fv_com)+exp(fv_reg))
untran_fv_reg <- exp(fv_reg)/(1+exp(fv_bjp)+exp(fv_com)+exp(fv_reg))
untran_fv_inc <- 1-(untran_fv_bjp + untran_fv_com + untran_fv_reg)

exp_bjp <- (untran_fv_bjp*4)-1
exp_com <- (untran_fv_com*4)-1
exp_reg <- (untran_fv_reg*4)-1
exp_inc <- (untran_fv_inc*4)-1

exp_bjpreg <- exp_bjp+exp_reg

hex_fvs <- terncoord(cbind(exp_bjpreg,exp_inc,exp_com))

#Full Hexagon plot - identify corners
l1 <- terncoords(c(1,-1,0))
l2 <- terncoords(c(-1,1,0))
l3 <- terncoords(c(1,0,-1))
l4 <- terncoords(c(-1,0,1))
l5 <- terncoords(c(0,-1,1))
l6 <- terncoords(c(0,1,-1))


##FIGURE 3
pdf(file= "Hex1.pdf", width = 5, height = 5, family = "Helvetica", pointsize = 10)
par(mfrow=c(1,1))
plot(hex_fvs[,1],hex_fvs[,2],xlim=c(-1.2,1.2),ylim=c(-1.2,1.2),pch=3,col="blue",cex=0.5,xlab="",ylab="",axes=FALSE)
points(x=l1[,1],y=l1[,2])
points(x=l2[,1],y=l2[,2])
points(x=l3[,1],y=l3[,2])
points(x=l4[,1],y=l4[,2])
points(x=l5[,1],y=l5[,2])
points(x=l6[,1],y=l6[,2])
polygon(x=c(-1.29,1.29,1.29,-1.29),y=c(1.29,1.29,-1.29,-1.29))


#Full hexagon contours - 20%iles
slope <- sqrt(3)
twentyunit <- sqrt(3)/5
for (i in -5:5){
  abline(i*twentyunit,slope,lty=2, col="darkgray")
  abline(i*twentyunit,-slope,lty=2, col="darkgray")
  abline(i*twentyunit/2,0,lty=2, col="darkgray")
}
abline(0*twentyunit,slope,lty=1,lwd=2)
abline(0*twentyunit,-slope,lty=1,lwd=2)
abline(0*twentyunit/2,0,lty=1,lwd=2)
lines(x=c(l1[,1],l2[,1]),y=c(l3[,2],l3[,2]),lty=1,lwd=2,col="green")
lines(x=c(l1[,1],l2[,1]),y=c(l4[,2],l4[,2]),lty=1,lwd=2,col="red")
lines(x=c(l1[,1],l5[,1]),y=c(l4[,2],l5[,2]),lty=1,lwd=2,col="green")
lines(x=c(l2[,1],l6[,1]),y=c(l1[,2],l5[,2]),lty=1,lwd=2,col="red")
lines(x=c(l6[,1],l3[,1]),y=c(l5[,2],l4[,2]),lty=1,lwd=2,col="green")
lines(x=c(l4[,1],l5[,1]),y=c(l3[,2],l5[,2]),lty=1,lwd=2,col="red")

#Interior hexagon contours
#Labels
text(0,1.1,"+ BJP/Regional")
text(-0.85,-0.5,"+ INC",srt=-57.5)
text(0.85,-0.5,"+ Far Left",srt=57.5)
text(0,-1.1,"- BJP/Regional")
text(+0.85,+0.5,"- INC",srt=-57.5)
text(-0.90,+0.5,"- Far Left",srt=57.5)
dev.off()


#############################################################

##Figure 4: Explaining the Hexagon

#Cool Plot
pdf(file= "Hexcool.pdf", width = 5, height = 5, family = "Helvetica", pointsize = 10)
plot(hex_fvs[,1],hex_fvs[,2],xlim=c(-1.2,1.2),ylim=c(-1.2,1.2),pch=3,col="white",cex=0.5,xlab="",ylab="",axes=FALSE)
points(x=l1[,1],y=l1[,2])
points(x=l2[,1],y=l2[,2])
points(x=l3[,1],y=l3[,2])
points(x=l4[,1],y=l4[,2])
points(x=l5[,1],y=l5[,2])
points(x=l6[,1],y=l6[,2])
polygon(x=c(-1.29,1.29,1.29,-1.29),y=c(1.29,1.29,-1.29,-1.29))

##Add Hexagon Contours
#Full hexagon contours - 20%iles
slope <- sqrt(3)
twentyunit <- sqrt(3)/5
for (i in -5:5){
  abline(i*twentyunit,slope,lty=2, col="gray")
  abline(i*twentyunit,-slope,lty=2, col="gray")
  abline(i*twentyunit/2,0,lty=2, col="gray")
}

##Add 0 Change Lines
abline(0*twentyunit,slope,lty=1,lwd=2)
abline(0*twentyunit,-slope,lty=1,lwd=2)
abline(0*twentyunit/2,0,lty=1,lwd=2)

##Add Hexagon Boundaries
lines(x=c(l1[,1],l2[,1]),y=c(l3[,2],l3[,2]),lty=1,lwd=2,col="green")
lines(x=c(l1[,1],l2[,1]),y=c(l4[,2],l4[,2]),lty=1,lwd=2,col="red")
lines(x=c(l1[,1],l5[,1]),y=c(l4[,2],l5[,2]),lty=1,lwd=2,col="green")
lines(x=c(l2[,1],l6[,1]),y=c(l1[,2],l5[,2]),lty=1,lwd=2,col="red")
lines(x=c(l6[,1],l3[,1]),y=c(l5[,2],l4[,2]),lty=1,lwd=2,col="green")
lines(x=c(l4[,1],l5[,1]),y=c(l3[,2],l5[,2]),lty=1,lwd=2,col="red")


#Label Hexagon Boundaries
par(font=2)
text(0,1,"+ BJP/Regional \n - INC \n - Far Left", cex=0.6)
text(-0.85,-0.5,"+ INC \n - BJP/Regional \n - Far Left",srt=-55.5, cex=0.6)
text(0.85,-0.5,"+ Far Left \n - INC \n - BJP/Regional",srt=55.5, cex=0.6)
text(0,-1.1,"- BJP/Regional \n + INC \n + Far Left", cex=0.6)
text(+0.85,+0.5,"- INC \n + BJP/Regional \n + Far Left",srt=-55.5, cex=0.6)
text(-0.90,+0.5,"- Far Left \n + INC \n + BJP/Regional",srt=55.5, cex=0.6)

##Add Labels
text(0,.6,"Pro- BJP", cex=1.2)
text(-0.5,-0.3,"Pro- INC", cex=1.2)
text(0.5,-0.3,"Pro- Far Left", cex=1.2)
text(0,-.6,"Pro- Left", cex=1.2)
text(+0.5,+0.25,"Polarized Voting", cex=1.2)
text(-0.5,+0.25,"Pro- Reform", cex=1.2)
dev.off()

#############################################################


##Hexagon Simplexes Figure 5: By Time Period

#####By time period #t=1 #Time FEs removed from models. 
m1_bjp_t1 <- lm(bjp_likeli2[data$t1==1] ~ tarmeasure[data$t1==1] + fdi_reform[data$t1==1] + licensing_reform[data$t1==1] + ps_reform[data$t1==1] + qr_reform[data$t1==1] + tardecreases[data$t1==1] + factor(region[data$t1==1]))
m1_com_t1 <- lm(com_likeli2[data$t1==1] ~ tarmeasure[data$t1==1] + fdi_reform[data$t1==1] + licensing_reform[data$t1==1] + ps_reform[data$t1==1] + qr_reform[data$t1==1] + tardecreases[data$t1==1] + factor(region[data$t1==1]))
m1_reg_t1 <- lm(reg_likeli2[data$t1==1] ~ tarmeasure[data$t1==1] + fdi_reform[data$t1==1] + licensing_reform[data$t1==1] + ps_reform[data$t1==1] + qr_reform[data$t1==1] + tardecreases[data$t1==1] + factor(region[data$t1==1]))

fv_bjp_t1 <- m1_bjp_t1$fitted.values
fv_com_t1 <- m1_com_t1$fitted.values
fv_reg_t1 <- m1_reg_t1$fitted.values
fv_inc_t1 <- 0 - (m1_bjp_t1$fitted.values + m1_com_t1$fitted.values + m1_reg_t1$fitted.values)

untran_fv_bjp_t1 <- exp(fv_bjp_t1)/(1+exp(fv_bjp_t1)+exp(fv_com_t1)+exp(fv_reg_t1))
untran_fv_com_t1 <- exp(fv_com_t1)/(1+exp(fv_bjp_t1)+exp(fv_com_t1)+exp(fv_reg_t1))
untran_fv_reg_t1 <- exp(fv_reg_t1)/(1+exp(fv_bjp_t1)+exp(fv_com_t1)+exp(fv_reg_t1))
untran_fv_inc_t1 <- 1-(untran_fv_bjp_t1 + untran_fv_com_t1 + untran_fv_reg_t1)

exp_bjp_t1 <- (untran_fv_bjp_t1*4)-1
exp_com_t1 <- (untran_fv_com_t1*4)-1
exp_reg_t1 <- (untran_fv_reg_t1*4)-1
exp_inc_t1 <- (untran_fv_inc_t1*4)-1

exp_bjpreg_t1 <- exp_bjp_t1+exp_reg_t1

hex_fvs_t1 <- terncoord(cbind(exp_bjpreg_t1,exp_inc_t1,exp_com_t1))

#Full Hexagon plot - identify corners
l1 <- terncoords(c(1,-1,0))
l2 <- terncoords(c(-1,1,0))
l3 <- terncoords(c(1,0,-1))
l4 <- terncoords(c(-1,0,1))
l5 <- terncoords(c(0,-1,1))
l6 <- terncoords(c(0,1,-1))


##FIGURE 5
pdf(file= "Hextime.pdf", width = 5, height = 5, family = "Helvetica", pointsize = 10)
par(mfrow=c(2,2),mar=c(1,1,1,1))
plot(hex_fvs_t1[,1],hex_fvs_t1[,2],xlim=c(-1.2,1.2),ylim=c(-1.2,1.2),pch=3,col="blue",cex=0.5,xlab="",ylab="",axes=FALSE)
points(x=l1[,1],y=l1[,2])
points(x=l2[,1],y=l2[,2])
points(x=l3[,1],y=l3[,2])
points(x=l4[,1],y=l4[,2])
points(x=l5[,1],y=l5[,2])
points(x=l6[,1],y=l6[,2])
polygon(x=c(-1.29,1.29,1.29,-1.29),y=c(1.29,1.29,-1.29,-1.29))

#Full hexagon contours - 20%iles
slope <- sqrt(3)
twentyunit <- sqrt(3)/5
for (i in -5:5){
  abline(i*twentyunit,slope,lty=2, col="darkgray")
  abline(i*twentyunit,-slope,lty=2, col="darkgray")
  abline(i*twentyunit/2,0,lty=2, col="darkgray")
}
abline(0*twentyunit,slope,lty=1,lwd=2)
abline(0*twentyunit,-slope,lty=1,lwd=2)
abline(0*twentyunit/2,0,lty=1,lwd=2)
lines(x=c(l1[,1],l2[,1]),y=c(l3[,2],l3[,2]),lty=1,lwd=2,col="green")
lines(x=c(l1[,1],l2[,1]),y=c(l4[,2],l4[,2]),lty=1,lwd=2,col="red")
lines(x=c(l1[,1],l5[,1]),y=c(l4[,2],l5[,2]),lty=1,lwd=2,col="green")
lines(x=c(l2[,1],l6[,1]),y=c(l1[,2],l5[,2]),lty=1,lwd=2,col="red")
lines(x=c(l6[,1],l3[,1]),y=c(l5[,2],l4[,2]),lty=1,lwd=2,col="green")
lines(x=c(l4[,1],l5[,1]),y=c(l3[,2],l5[,2]),lty=1,lwd=2,col="red")

#Interior hexagon contours
#Labels
text(0,1.1,"+ BJP/Regional")
text(-0.85,-0.5,"+ INC",srt=-57.5)
text(0.85,-0.5,"+ Far Left",srt=57.5)
text(0,-1.1,"- BJP/Regional")
text(+0.85,+0.5,"- INC",srt=-57.5)
text(-0.90,+0.5,"- Far Left",srt=57.5)
text(-1,+1.1,"1991-96")


#####By time period #t=2#Time FEs removed from models. 
m1_bjp_t2 <- lm(bjp_likeli2[data$t2==1] ~ tarmeasure[data$t2==1] + fdi_reform[data$t2==1] + licensing_reform[data$t2==1] + ps_reform[data$t2==1] + qr_reform[data$t2==1] + tardecreases[data$t2==1] + factor(region[data$t2==1]))
m1_com_t2 <- lm(com_likeli2[data$t2==1] ~ tarmeasure[data$t2==1] + fdi_reform[data$t2==1] + licensing_reform[data$t2==1] + ps_reform[data$t2==1] + qr_reform[data$t2==1] + tardecreases[data$t2==1] + factor(region[data$t2==1]))
m1_reg_t2 <- lm(reg_likeli2[data$t2==1] ~ tarmeasure[data$t2==1] + fdi_reform[data$t2==1] + licensing_reform[data$t2==1] + ps_reform[data$t2==1] + qr_reform[data$t2==1] + tardecreases[data$t2==1] + factor(region[data$t2==1]))

fv_bjp_t2 <- m1_bjp_t2$fitted.values
fv_com_t2 <- m1_com_t2$fitted.values
fv_reg_t2 <- m1_reg_t2$fitted.values
fv_inc_t2 <- 0 - (m1_bjp_t2$fitted.values + m1_com_t2$fitted.values + m1_reg_t2$fitted.values)

untran_fv_bjp_t2 <- exp(fv_bjp_t2)/(1+exp(fv_bjp_t2)+exp(fv_com_t2)+exp(fv_reg_t2))
untran_fv_com_t2 <- exp(fv_com_t2)/(1+exp(fv_bjp_t2)+exp(fv_com_t2)+exp(fv_reg_t2))
untran_fv_reg_t2 <- exp(fv_reg_t2)/(1+exp(fv_bjp_t2)+exp(fv_com_t2)+exp(fv_reg_t2))
untran_fv_inc_t2 <- 1-(untran_fv_bjp_t2 + untran_fv_com_t2 + untran_fv_reg_t2)

exp_bjp_t2 <- (untran_fv_bjp_t2*4)-1
exp_com_t2 <- (untran_fv_com_t2*4)-1
exp_reg_t2 <- (untran_fv_reg_t2*4)-1
exp_inc_t2 <- (untran_fv_inc_t2*4)-1

exp_bjpreg_t2 <- exp_bjp_t2+exp_reg_t2

hex_fvs_t2 <- terncoord(cbind(exp_bjpreg_t2,exp_inc_t2,exp_com_t2))

#Full Hexagon plot - identify corners
l1 <- terncoords(c(1,-1,0))
l2 <- terncoords(c(-1,1,0))
l3 <- terncoords(c(1,0,-1))
l4 <- terncoords(c(-1,0,1))
l5 <- terncoords(c(0,-1,1))
l6 <- terncoords(c(0,1,-1))

plot(hex_fvs_t2[,1],hex_fvs_t2[,2],xlim=c(-1.2,1.2),ylim=c(-1.2,1.2),pch=3,col="blue",cex=0.5,xlab="",ylab="",axes=FALSE)
points(x=l1[,1],y=l1[,2])
points(x=l2[,1],y=l2[,2])
points(x=l3[,1],y=l3[,2])
points(x=l4[,1],y=l4[,2])
points(x=l5[,1],y=l5[,2])
points(x=l6[,1],y=l6[,2])
polygon(x=c(-1.29,1.29,1.29,-1.29),y=c(1.29,1.29,-1.29,-1.29))

#Full hexagon contours - 20%iles
slope <- sqrt(3)
twentyunit <- sqrt(3)/5
for (i in -5:5){
  abline(i*twentyunit,slope,lty=2, col="darkgray")
  abline(i*twentyunit,-slope,lty=2, col="darkgray")
  abline(i*twentyunit/2,0,lty=2, col="darkgray")
}
abline(0*twentyunit,slope,lty=1,lwd=2)
abline(0*twentyunit,-slope,lty=1,lwd=2)
abline(0*twentyunit/2,0,lty=1,lwd=2)
lines(x=c(l1[,1],l2[,1]),y=c(l3[,2],l3[,2]),lty=1,lwd=2,col="green")
lines(x=c(l1[,1],l2[,1]),y=c(l4[,2],l4[,2]),lty=1,lwd=2,col="red")
lines(x=c(l1[,1],l5[,1]),y=c(l4[,2],l5[,2]),lty=1,lwd=2,col="green")
lines(x=c(l2[,1],l6[,1]),y=c(l1[,2],l5[,2]),lty=1,lwd=2,col="red")
lines(x=c(l6[,1],l3[,1]),y=c(l5[,2],l4[,2]),lty=1,lwd=2,col="green")
lines(x=c(l4[,1],l5[,1]),y=c(l3[,2],l5[,2]),lty=1,lwd=2,col="red")

#Interior hexagon contours
#Labels
text(0,1.1,"+ BJP/Regional")
text(-0.85,-0.5,"+ INC",srt=-57.5)
text(0.85,-0.5,"+ Far Left",srt=57.5)
text(0,-1.1,"- BJP/Regional")
text(+0.85,+0.5,"- INC",srt=-57.5)
text(-0.90,+0.5,"- Far Left",srt=57.5)
text(-1,+1.1,"1996-98")

#####By time period #t=3#Time FEs removed from models.

m1_bjp_t3 <- lm(bjp_likeli2[data$t3==1] ~ tarmeasure[data$t3==1] + fdi_reform[data$t3==1] + licensing_reform[data$t3==1] + ps_reform[data$t3==1] + qr_reform[data$t3==1] + tardecreases[data$t3==1] + factor(region[data$t3==1]))
m1_com_t3 <- lm(com_likeli2[data$t3==1] ~ tarmeasure[data$t3==1] + fdi_reform[data$t3==1] + licensing_reform[data$t3==1] + ps_reform[data$t3==1] + qr_reform[data$t3==1] + tardecreases[data$t3==1] + factor(region[data$t3==1]))
m1_reg_t3 <- lm(reg_likeli2[data$t3==1] ~ tarmeasure[data$t3==1] + fdi_reform[data$t3==1] + licensing_reform[data$t3==1] + ps_reform[data$t3==1] + qr_reform[data$t3==1] + tardecreases[data$t3==1] + factor(region[data$t3==1]))

fv_bjp_t3 <- m1_bjp_t3$fitted.values
fv_com_t3 <- m1_com_t3$fitted.values
fv_reg_t3 <- m1_reg_t3$fitted.values
fv_inc_t3 <- 0 - (m1_bjp_t3$fitted.values + m1_com_t3$fitted.values + m1_reg_t3$fitted.values)

untran_fv_bjp_t3 <- exp(fv_bjp_t3)/(1+exp(fv_bjp_t3)+exp(fv_com_t3)+exp(fv_reg_t3))
untran_fv_com_t3 <- exp(fv_com_t3)/(1+exp(fv_bjp_t3)+exp(fv_com_t3)+exp(fv_reg_t3))
untran_fv_reg_t3 <- exp(fv_reg_t3)/(1+exp(fv_bjp_t3)+exp(fv_com_t3)+exp(fv_reg_t3))
untran_fv_inc_t3 <- 1-(untran_fv_bjp_t3 + untran_fv_com_t3 + untran_fv_reg_t3)

exp_bjp_t3 <- (untran_fv_bjp_t3*4)-1
exp_com_t3 <- (untran_fv_com_t3*4)-1
exp_reg_t3 <- (untran_fv_reg_t3*4)-1
exp_inc_t3 <- (untran_fv_inc_t3*4)-1

exp_bjpreg_t3 <- exp_bjp_t3+exp_reg_t3

hex_fvs_t3 <- terncoord(cbind(exp_bjpreg_t3,exp_inc_t3,exp_com_t3))

#Full Hexagon plot - identify corners
l1 <- terncoords(c(1,-1,0))
l2 <- terncoords(c(-1,1,0))
l3 <- terncoords(c(1,0,-1))
l4 <- terncoords(c(-1,0,1))
l5 <- terncoords(c(0,-1,1))
l6 <- terncoords(c(0,1,-1))

plot(hex_fvs_t3[,1],hex_fvs_t3[,2],xlim=c(-1.2,1.2),ylim=c(-1.2,1.2),pch=3,col="blue",cex=0.5,xlab="",ylab="",axes=FALSE)
points(x=l1[,1],y=l1[,2])
points(x=l2[,1],y=l2[,2])
points(x=l3[,1],y=l3[,2])
points(x=l4[,1],y=l4[,2])
points(x=l5[,1],y=l5[,2])
points(x=l6[,1],y=l6[,2])
polygon(x=c(-1.29,1.29,1.29,-1.29),y=c(1.29,1.29,-1.29,-1.29))

#Full hexagon contours - 20%iles
slope <- sqrt(3)
twentyunit <- sqrt(3)/5
for (i in -5:5){
  abline(i*twentyunit,slope,lty=2, col="darkgray")
  abline(i*twentyunit,-slope,lty=2, col="darkgray")
  abline(i*twentyunit/2,0,lty=2, col="darkgray")
}
abline(0*twentyunit,slope,lty=1,lwd=2)
abline(0*twentyunit,-slope,lty=1,lwd=2)
abline(0*twentyunit/2,0,lty=1,lwd=2)
lines(x=c(l1[,1],l2[,1]),y=c(l3[,2],l3[,2]),lty=1,lwd=2,col="green")
lines(x=c(l1[,1],l2[,1]),y=c(l4[,2],l4[,2]),lty=1,lwd=2,col="red")
lines(x=c(l1[,1],l5[,1]),y=c(l4[,2],l5[,2]),lty=1,lwd=2,col="green")
lines(x=c(l2[,1],l6[,1]),y=c(l1[,2],l5[,2]),lty=1,lwd=2,col="red")
lines(x=c(l6[,1],l3[,1]),y=c(l5[,2],l4[,2]),lty=1,lwd=2,col="green")
lines(x=c(l4[,1],l5[,1]),y=c(l3[,2],l5[,2]),lty=1,lwd=2,col="red")

#Interior hexagon contours
#Labels
text(0,1.1,"+ BJP/Regional")
text(-0.85,-0.5,"+ INC",srt=-57.5)
text(0.85,-0.5,"+ Far Left",srt=57.5)
text(0,-1.1,"- BJP/Regional")
text(+0.85,+0.5,"- INC",srt=-57.5)
text(-0.90,+0.5,"- Far Left",srt=57.5)
text(-1,+1.1,"1998-99")

#####By time period #t=4#Time FEs removed from models.
m1_bjp_t4 <- lm(bjp_likeli2[data$t4==1] ~ tarmeasure[data$t4==1] + fdi_reform[data$t4==1] + licensing_reform[data$t4==1] + ps_reform[data$t4==1] + qr_reform[data$t4==1] + tardecreases[data$t4==1] + factor(region[data$t4==1]))
m1_com_t4 <- lm(com_likeli2[data$t4==1] ~ tarmeasure[data$t4==1] + fdi_reform[data$t4==1] + licensing_reform[data$t4==1] + ps_reform[data$t4==1] + qr_reform[data$t4==1] + tardecreases[data$t4==1] + factor(region[data$t4==1]))
m1_reg_t4 <- lm(reg_likeli2[data$t4==1] ~ tarmeasure[data$t4==1] + fdi_reform[data$t4==1] + licensing_reform[data$t4==1] + ps_reform[data$t4==1] + qr_reform[data$t4==1] + tardecreases[data$t4==1] + factor(region[data$t4==1]))

fv_bjp_t4 <- m1_bjp_t4$fitted.values
fv_com_t4 <- m1_com_t4$fitted.values
fv_reg_t4 <- m1_reg_t4$fitted.values
fv_inc_t4 <- 0 - (m1_bjp_t4$fitted.values + m1_com_t4$fitted.values + m1_reg_t4$fitted.values)

untran_fv_bjp_t4 <- exp(fv_bjp_t4)/(1+exp(fv_bjp_t4)+exp(fv_com_t4)+exp(fv_reg_t4))
untran_fv_com_t4 <- exp(fv_com_t4)/(1+exp(fv_bjp_t4)+exp(fv_com_t4)+exp(fv_reg_t4))
untran_fv_reg_t4 <- exp(fv_reg_t4)/(1+exp(fv_bjp_t4)+exp(fv_com_t4)+exp(fv_reg_t4))
untran_fv_inc_t4 <- 1-(untran_fv_bjp_t4 + untran_fv_com_t4 + untran_fv_reg_t4)

exp_bjp_t4 <- (untran_fv_bjp_t4*4)-1
exp_com_t4 <- (untran_fv_com_t4*4)-1
exp_reg_t4 <- (untran_fv_reg_t4*4)-1
exp_inc_t4 <- (untran_fv_inc_t4*4)-1

exp_bjpreg_t4 <- exp_bjp_t4+exp_reg_t4

hex_fvs_t4 <- terncoord(cbind(exp_bjpreg_t4,exp_inc_t4,exp_com_t4))

#Full Hexagon plot - identify corners
l1 <- terncoords(c(1,-1,0))
l2 <- terncoords(c(-1,1,0))
l3 <- terncoords(c(1,0,-1))
l4 <- terncoords(c(-1,0,1))
l5 <- terncoords(c(0,-1,1))
l6 <- terncoords(c(0,1,-1))

plot(hex_fvs_t4[,1],hex_fvs_t4[,2],xlim=c(-1.2,1.2),ylim=c(-1.2,1.2),pch=3,col="blue",cex=0.5,xlab="",ylab="",axes=FALSE)
points(x=l1[,1],y=l1[,2])
points(x=l2[,1],y=l2[,2])
points(x=l3[,1],y=l3[,2])
points(x=l4[,1],y=l4[,2])
points(x=l5[,1],y=l5[,2])
points(x=l6[,1],y=l6[,2])
polygon(x=c(-1.29,1.29,1.29,-1.29),y=c(1.29,1.29,-1.29,-1.29))

#Full hexagon contours - 20%iles
slope <- sqrt(3)
twentyunit <- sqrt(3)/5
for (i in -5:5){
  abline(i*twentyunit,slope,lty=2, col="darkgray")
  abline(i*twentyunit,-slope,lty=2, col="darkgray")
  abline(i*twentyunit/2,0,lty=2, col="darkgray")
}
abline(0*twentyunit,slope,lty=1,lwd=2)
abline(0*twentyunit,-slope,lty=1,lwd=2)
abline(0*twentyunit/2,0,lty=1,lwd=2)
lines(x=c(l1[,1],l2[,1]),y=c(l3[,2],l3[,2]),lty=1,lwd=2,col="green")
lines(x=c(l1[,1],l2[,1]),y=c(l4[,2],l4[,2]),lty=1,lwd=2,col="red")
lines(x=c(l1[,1],l5[,1]),y=c(l4[,2],l5[,2]),lty=1,lwd=2,col="green")
lines(x=c(l2[,1],l6[,1]),y=c(l1[,2],l5[,2]),lty=1,lwd=2,col="red")
lines(x=c(l6[,1],l3[,1]),y=c(l5[,2],l4[,2]),lty=1,lwd=2,col="green")
lines(x=c(l4[,1],l5[,1]),y=c(l3[,2],l5[,2]),lty=1,lwd=2,col="red")

#Interior hexagon contours
#Labels
text(0,1.1,"+ BJP/Regional")
text(-0.85,-0.5,"+ INC",srt=-57.5)
text(0.85,-0.5,"+ Far Left",srt=57.5)
text(0,-1.1,"- BJP/Regional")
text(+0.85,+0.5,"- INC",srt=-57.5)
text(-0.90,+0.5,"- Far Left",srt=57.5)
text(-1,+1.1,"1999-04")
dev.off()


#############################################################

##FIGURE NOT INCLUDED: First Differences

bjpreg_fd <- fd2_bjp + fd2_reg
inc_fd <- fd2_inc
com_fd <- fd2_inc
fd <- cbind(bjpreg_fd,inc_fd,com_fd)

hex_fd <- terncoord(fd)

par(mfrow=c(1,1))
pdf(file= "Hexfd.pdf", width = 5, height = 5, family = "Helvetica", pointsize = 10)
plot(hex_fd[,1],hex_fd[,2],xlim=c(-1.2,1.2),ylim=c(-1.2,1.2),pch=3,col="blue",cex=0.5,xlab="",ylab="",axes=FALSE)
points(x=l1[,1],y=l1[,2])
points(x=l2[,1],y=l2[,2])
points(x=l3[,1],y=l3[,2])
points(x=l4[,1],y=l4[,2])
points(x=l5[,1],y=l5[,2])
points(x=l6[,1],y=l6[,2])
polygon(x=c(-1.29,1.29,1.29,-1.29),y=c(1.29,1.29,-1.29,-1.29))


#Full hexagon contours - 20%iles
slope <- sqrt(3)
twentyunit <- sqrt(3)/5
for (i in -5:5){
  abline(i*twentyunit,slope,lty=2, col="darkgray")
  abline(i*twentyunit,-slope,lty=2, col="darkgray")
  abline(i*twentyunit/2,0,lty=2, col="darkgray")
}
abline(0*twentyunit,slope,lty=1,lwd=2)
abline(0*twentyunit,-slope,lty=1,lwd=2)
abline(0*twentyunit/2,0,lty=1,lwd=2)
lines(x=c(l1[,1],l2[,1]),y=c(l3[,2],l3[,2]),lty=1,lwd=2,col="green")
lines(x=c(l1[,1],l2[,1]),y=c(l4[,2],l4[,2]),lty=1,lwd=2,col="red")
lines(x=c(l1[,1],l5[,1]),y=c(l4[,2],l5[,2]),lty=1,lwd=2,col="green")
lines(x=c(l2[,1],l6[,1]),y=c(l1[,2],l5[,2]),lty=1,lwd=2,col="red")
lines(x=c(l6[,1],l3[,1]),y=c(l5[,2],l4[,2]),lty=1,lwd=2,col="green")
lines(x=c(l4[,1],l5[,1]),y=c(l3[,2],l5[,2]),lty=1,lwd=2,col="red")

#Interior hexagon contours
#Labels
text(0,1.1,"+ BJP/Regional")
text(-0.85,-0.5,"+ INC",srt=-57.5)
text(0.85,-0.5,"+ Far Left",srt=57.5)
text(0,-1.1,"- BJP/Regional")
text(+0.85,+0.5,"- INC",srt=-57.5)
text(-0.90,+0.5,"- Far Left",srt=57.5)
dev.off()


#############################################################
##----------------------Section 5--------------------------##
##--------------Reanalyzing Tandon (2012)------------------##
#############################################################

##Generate Lags
data.plm$incum_1<-lag(data.plm$incum_voteshare,1)
data.plm$left_1<-lag(data.plm$left,1)
data.plm$bjp_1<-lag(data.plm$bjp,1)
data.plm$inc_1<-lag(data.plm$inc,1)
data.plm$com_1<-lag(data.plm$com,1)
data.plm$regional_1<-lag(data.plm$regional,1)

lagt1 <- seq(from=1,to=472,by=4)
lagt2 <- seq(from=1,to=472,by=4)+1
lagt3 <- seq(from=1,to=472,by=4)+2
lagt4 <- seq(from=1,to=472,by=4)+3

bjp.t1<-data.plm$bjp_1[lagt1]
bjp.t2<-data.plm$bjp_1[lagt2]
bjp.t3<-data.plm$bjp_1[lagt3]
bjp.t4<-data.plm$bjp_1[lagt4]
bjp.lag<-c(bjp.t1,bjp.t2, bjp.t3, bjp.t4)

com.t1<-data.plm$com_1[lagt1]
com.t2<-data.plm$com_1[lagt2]
com.t3<-data.plm$com_1[lagt3]
com.t4<-data.plm$com_1[lagt4]
com.lag<-c(com.t1,com.t2, com.t3, com.t4)

reg.t1<-data.plm$regional_1[lagt1]
reg.t2<-data.plm$regional_1[lagt2]
reg.t3<-data.plm$regional_1[lagt3]
reg.t4<-data.plm$regional_1[lagt4]
reg.lag<-c(reg.t1,reg.t2, reg.t3, reg.t4)


#Run separate regressions using OLS
m1_bjp_2 <- lm(bjp_likeli2 ~ tarmeasure + fdi_reform + licensing_reform + ps_reform + qr_reform + tardecreases + bjp.lag + factor(time) + factor(region))
m1_com_2 <- lm(com_likeli2 ~ tarmeasure + fdi_reform + licensing_reform + ps_reform + qr_reform + tardecreases + com.lag + factor(time) + factor(region))
m1_reg_2 <- lm(reg_likeli2 ~ tarmeasure + fdi_reform + licensing_reform + ps_reform + qr_reform + tardecreases + reg.lag + factor(time) + factor(region))
summary(m1_bjp_2)

#############################################################

#SUR model

##Subset Data
bjp.like<-bjp_likeli2[119:472]
com.like<-com_likeli2[119:472]
reg.like<-reg_likeli2[119:472]

comp.sub<-comp[t2==1 | t3==1 | t4==1,]
lags<-na.omit(cbind(bjp.lag, com.lag, reg.lag))


library(systemfit)
one_2 <- bjp.like ~ comp.sub$tarmeasure + comp.sub$fdi_reform + comp.sub$licensing_reform + comp.sub$ps_reform + comp.sub$tardecreases + lags[,1] + factor(comp.sub$time) + factor(comp.sub$region)
two_2 <- com.like ~ comp.sub$tarmeasure + comp.sub$fdi_reform + comp.sub$licensing_reform + comp.sub$ps_reform + comp.sub$tardecreases + lags[,2] + factor(comp.sub$time) + factor(comp.sub$region)
three_2 <- reg.like ~ comp.sub$tarmeasure + comp.sub$fdi_reform + comp.sub$licensing_reform + comp.sub$ps_reform  + comp.sub$tardecreases + lags[,3] + factor(comp.sub$time) + factor(comp.sub$region)
form_2 <- list(one_2, two_2, three_2)
SURM_2 <- systemfit(form_2,"SUR")
summary(SURM_2)
names(SURM_2)
vcov.SURM_2<-vcov(SURM_2)


#############################################################

##1. Expected Values

##Set All Covariates at their median
median_bjp<-c(1,apply(cbind(comp.sub[,13],comp.sub[,16:18], comp.sub[,15], lags[,1], comp.sub[,5:6]),2,FUN=median), rep(0,117))
median_com<-c(1,apply(cbind(comp.sub[,13],comp.sub[,16:18], comp.sub[,15], lags[,2], comp.sub[,5:6]),2,FUN=median), rep(0,117))
median_reg<-c(1,apply(cbind(comp.sub[,13],comp.sub[,16:18], comp.sub[,15], lags[,3], comp.sub[,5:6]),2,FUN=median), rep(0,117))

##BJP
##Simulate Betas for BJP Model
set.seed(1234)
simbetasm1_bjp_2 <- rmvnorm(1000,mean=SURM_2$coefficients[1:126],sigma=vcov.SURM_2[1:126,1:126])

##Estimate Expected Values
ev.m1_bjp_2<-simbetasm1_bjp_2 %*% median_bjp
summary(ev.m1_bjp_2)


##COM
##Simulate Betas for COM Model
set.seed(1234)
simbetasm1_com_2 <- rmvnorm(1000,mean=SURM_2$coefficients[127:252],sigma=vcov.SURM_2[127:252,127:252])

##Estimate Expected Values
ev.m1_com_2<-simbetasm1_com_2 %*% median_com
summary(ev.m1_com_2)


##REGIONAL
##Simulate Betas for Regional Model
set.seed(1234)
simbetasm1_reg_2 <- rmvnorm(1000,mean=SURM_2$coefficients[253:378],sigma=vcov.SURM_2[253:378,253:378])

##Estimate Expected Values
ev.m1_reg_2<-simbetasm1_reg_2 %*% median_reg
summary(ev.m1_reg_2)


#Untransform
untran_bjp_ldv <- exp(ev.m1_bjp_2)/(1+exp(ev.m1_bjp_2)+exp(ev.m1_com_2)+exp(ev.m1_reg_2))
untran_com_ldv <- exp(ev.m1_com_2)/(1+exp(ev.m1_bjp_2)+exp(ev.m1_com_2)+exp(ev.m1_reg_2))
untran_reg_ldv <- exp(ev.m1_reg_2)/(1+exp(ev.m1_bjp_2)+exp(ev.m1_com_2)+exp(ev.m1_reg_2))
untran_inc_ldv <- 1-(untran_bjp_ldv + untran_com_ldv + untran_reg_ldv)

exp_bjp_ldv <- (untran_bjp_ldv*4)-1
exp_com_ldv <- (untran_com_ldv*4)-1
exp_reg_ldv <- (untran_reg_ldv*4)-1
exp_inc_ldv <- (untran_inc_ldv*4)-1

expected.values_ldv<-cbind(exp_inc_ldv, exp_bjp_ldv, exp_com_ldv, exp_reg_ldv)
colnames(expected.values_ldv)<-c("EV[INC]", "EV[BJP]", "EV[COM]", "EV[REG]")
summary(expected.values_ldv)



##FIGURE 6: Plot Expected Values
pdf("ev2.pdf")
par(mfrow=c(1,1))
plot(density(expected.values_ldv[,3]), col="gold2", xlim=c(-.5, .5), lwd=2, xlab="Simulated Expected Values", sub="Covariates at Medians", main="Expected Values using Compositional \n Voting Models and an LDV")
lines(density(expected.values_ldv[,1]), col="magenta4", lwd=2)
lines(density(expected.values_ldv[,2]), col="darkseagreen4", lwd=2)
lines(density(expected.values_ldv[,4]), col="orangered3", lwd=2)
par(font=2)
legend(x="topright", legend=c("BJP","INC", "Far Left", "Regional"), col=c("darkseagreen4", "magenta4", "gold2", "orangered3"), lty=c(1,1,1,1), lwd=c(2,2,2,2))
dev.off()

#############################################################

##2. FIRST DIFFERENCES

#Set values of Tariff Measue at 0 and minimum
notar_bjp <- median_bjp
notar_bjp[2] <- 0
tarmin_bjp <- median_bjp
tarmin_bjp[2] <- min(tarmeasure)

notar_com <- median_com
notar_com[2] <- 0
tarmin_com <- median_com
tarmin_com[2] <- min(tarmeasure)

notar_reg <- median_reg
notar_reg[2] <- 0
tarmin_reg <- median_reg
tarmin_reg[2] <- min(tarmeasure)

#RCalculate Expected Values Conditional on Both Levels of Tariffs
notarbjp_2 <- (simbetasm1_bjp_2 %*% notar_bjp)
tarbjp_2 <- (simbetasm1_bjp_2 %*% tarmin_bjp)
notarcom_2 <- (simbetasm1_com_2 %*% notar_com)
tarcom_2 <- (simbetasm1_com_2 %*% tarmin_com)
notarreg_2 <- (simbetasm1_reg_2 %*% notar_reg)
tarreg_2 <- (simbetasm1_reg_2 %*% tarmin_reg)

##Untransform Expected Values
untran_notarcom_2 <- exp(notarcom_2)/(1+exp(notarbjp_2)+exp(notarcom_2)+exp(notarreg_2))
untran_notarreg_2 <- exp(notarreg_2)/(1+exp(notarbjp_2)+exp(notarcom_2)+exp(notarreg_2))
untran_notarbjp_2 <- exp(notarbjp_2)/(1+exp(notarbjp_2)+exp(notarcom_2)+exp(notarreg_2))
untran_notarinc_2 <- 1-(untran_notarbjp_2 + untran_notarcom_2 + untran_notarreg_2)

notar_fd_bjp_2 <- (untran_notarbjp_2*4)-1
notar_fd_com_2 <- (untran_notarcom_2*4)-1
notar_fd_reg_2 <- (untran_notarreg_2*4)-1
notar_fd_inc_2 <- (untran_notarinc_2*4)-1

untran_tarcom_2 <- exp(tarcom_2)/(1+exp(tarbjp_2)+exp(tarcom_2)+exp(tarreg_2))
untran_tarreg_2 <- exp(tarreg_2)/(1+exp(tarbjp_2)+exp(tarcom_2)+exp(tarreg_2))
untran_tarbjp_2 <- exp(tarbjp_2)/(1+exp(tarbjp_2)+exp(tarcom_2)+exp(tarreg_2))
untran_tarinc_2 <- 1-(untran_tarbjp_2 + untran_tarcom_2 + untran_tarreg_2)

tar_fd_bjp_2 <- (untran_tarbjp_2*4)-1
tar_fd_com_2 <- (untran_tarcom_2*4)-1
tar_fd_reg_2 <- (untran_tarreg_2*4)-1
tar_fd_inc_2 <- (untran_tarinc_2*4)-1


##Calculate First Differences Given Untransformed Values
fd2_bjp_2 <- notar_fd_bjp_2 - tar_fd_bjp_2
fd2_com_2 <- notar_fd_com_2 - tar_fd_com_2
fd2_reg_2 <- notar_fd_reg_2 - tar_fd_reg_2
fd2_inc_2 <- notar_fd_inc_2 - tar_fd_inc_2

first.differences_ldv<-cbind(fd2_inc_2, fd2_bjp_2, fd2_com_2, fd2_reg_2)
colnames(first.differences_ldv)<-c("FD[INC]", "FD[BJP]", "FD[COM]", "FD[REG]")
summary(first.differences_ldv)



##FIGURE 7: Plot First Differences
pdf("fd2.pdf")
plot(density(first.differences_ldv[,3]), col="gold2", xlim=c(-0.7, 0.7), lwd=2, xlab="Simulated First Difference", sub="EV[Vote Change | No Tariff] - EV[Vote Change | Min Tariff]", main="First Differences using Compositional \n Voting Models and an LDV")
lines(density(first.differences_ldv[,1]), col="magenta4", lwd=2)
lines(density(first.differences_ldv[,2]), col="darkseagreen4", lwd=2)
lines(density(first.differences_ldv[,4]), col="orangered3", lwd=2)
par(font=2)
legend(x="topright", legend=c("BJP","INC", "Far Left", "Regional"), col=c("darkseagreen4", "magenta4", "gold2", "orangered3"), lty=c(1,1,1,1), lwd=c(2,2,2,2))
dev.off()



#############################################################
##-------------Plotting Hexagonal Simplexes----------------##
#############################################################

#Predicted values in original scale
fv_bjp_ldv <- m1_bjp_2$fitted.values
fv_com_ldv <- m1_com_2$fitted.values
fv_reg_ldv <- m1_reg_2$fitted.values
fv_inc_ldv <- 0 - (m1_bjp_2$fitted.values + m1_com_2$fitted.values + m1_reg_2$fitted.values)

untran_fv_bjp_ldv <- exp(fv_bjp_ldv)/(1+exp(fv_bjp_ldv)+exp(fv_com_ldv)+exp(fv_reg_ldv))
untran_fv_com_ldv <- exp(fv_com_ldv)/(1+exp(fv_bjp_ldv)+exp(fv_com_ldv)+exp(fv_reg_ldv))
untran_fv_reg_ldv <- exp(fv_reg_ldv)/(1+exp(fv_bjp_ldv)+exp(fv_com_ldv)+exp(fv_reg_ldv))
untran_fv_inc_ldv <- 1-(untran_fv_bjp_ldv + untran_fv_com_ldv + untran_fv_reg_ldv)

exp_bjp_ldv <- (untran_fv_bjp_ldv*4)-1
exp_com_ldv <- (untran_fv_com_ldv*4)-1
exp_reg_ldv <- (untran_fv_reg_ldv*4)-1
exp_inc_ldv <- (untran_fv_inc_ldv*4)-1

exp_bjpreg_ldv <- exp_bjp_ldv+exp_reg_ldv

hex_fvs_ldv <- terncoord(cbind(exp_bjpreg_ldv,exp_inc_ldv,exp_com_ldv))


##FIGURE NOT INCLUDED: Fitted Values
pdf(file= "Hex1_ldv.pdf", width = 5, height = 5, family = "Helvetica", pointsize = 10)
plot(hex_fvs_ldv[,1],hex_fvs_ldv[,2],xlim=c(-1.2,1.2),ylim=c(-1.2,1.2),pch=3,col="blue",cex=0.5,xlab="",ylab="",axes=FALSE)
points(x=l1[,1],y=l1[,2])
points(x=l2[,1],y=l2[,2])
points(x=l3[,1],y=l3[,2])
points(x=l4[,1],y=l4[,2])
points(x=l5[,1],y=l5[,2])
points(x=l6[,1],y=l6[,2])
polygon(x=c(-1.29,1.29,1.29,-1.29),y=c(1.29,1.29,-1.29,-1.29))

#Full hexagon contours - 20%iles
slope <- sqrt(3)
twentyunit <- sqrt(3)/5
for (i in -5:5){
  abline(i*twentyunit,slope,lty=2, col="darkgray")
  abline(i*twentyunit,-slope,lty=2, col="darkgray")
  abline(i*twentyunit/2,0,lty=2, col="darkgray")
}
abline(0*twentyunit,slope,lty=1,lwd=2)
abline(0*twentyunit,-slope,lty=1,lwd=2)
abline(0*twentyunit/2,0,lty=1,lwd=2)
lines(x=c(l1[,1],l2[,1]),y=c(l3[,2],l3[,2]),lty=1,lwd=2,col="green")
lines(x=c(l1[,1],l2[,1]),y=c(l4[,2],l4[,2]),lty=1,lwd=2,col="red")
lines(x=c(l1[,1],l5[,1]),y=c(l4[,2],l5[,2]),lty=1,lwd=2,col="green")
lines(x=c(l2[,1],l6[,1]),y=c(l1[,2],l5[,2]),lty=1,lwd=2,col="red")
lines(x=c(l6[,1],l3[,1]),y=c(l5[,2],l4[,2]),lty=1,lwd=2,col="green")
lines(x=c(l4[,1],l5[,1]),y=c(l3[,2],l5[,2]),lty=1,lwd=2,col="red")

#Interior hexagon contours
#Labels
text(0,1.1,"+ BJP/Regional")
text(-0.85,-0.5,"+ INC",srt=-57.5)
text(0.85,-0.5,"+ Far Left",srt=57.5)
text(0,-1.1,"- BJP/Regional")
text(+0.85,+0.5,"- INC",srt=-57.5)
text(-0.90,+0.5,"- Far Left",srt=57.5)
dev.off()

#############################################################

##FIGURE NOT INCLUDED: By Time Period Fitted Values

#####By time period #t=1 #Time FEs removed from models. 
lagt1 <- seq(from=1,to=472,by=4)
lagt2 <- seq(from=1,to=472,by=4)+1
lagt3 <- seq(from=1,to=472,by=4)+2
lagt4 <- seq(from=1,to=472,by=4)+3

#Correctly, t1 fails as lagt1 is null
#m1_bjp_2_t1 <- lm(bjp_likeli2[data$t1==1] ~ tarmeasure[data$t1==1] + fdi_reform[data$t1==1] + licensing_reform[data$t1==1] + ps_reform[data$t1==1] + qr_reform[data$t1==1] + tardecreases[data$t1==1] + bjp.t1 + factor(region[data$t1==1]))
#m1_com_2_t1 <- lm(com_likeli2[data$t1==1] ~ tarmeasure[data$t1==1] + fdi_reform[data$t1==1] + licensing_reform[data$t1==1] + ps_reform[data$t1==1] + qr_reform[data$t1==1] + tardecreases[data$t1==1] + com.t1 + factor(region[data$t1==1]))
#m1_reg_2_t1 <- lm(reg_likeli2[data$t1==1] ~ tarmeasure[data$t1==1] + fdi_reform[data$t1==1] + licensing_reform[data$t1==1] + ps_reform[data$t1==1] + qr_reform[data$t1==1] + tardecreases[data$t1==1] + reg.t1 + factor(region[data$t1==1]))

#####By time period #t=2 #Time FEs removed from models.
m1_bjp_2_t2 <- lm(bjp_likeli2[data$t2==1] ~ tarmeasure[data$t2==1] + fdi_reform[data$t2==1] + licensing_reform[data$t2==1] + ps_reform[data$t2==1] + qr_reform[data$t2==1] + tardecreases[data$t2==1] + bjp.t2 + factor(region[data$t2==1]))
m1_com_2_t2 <- lm(com_likeli2[data$t2==1] ~ tarmeasure[data$t2==1] + fdi_reform[data$t2==1] + licensing_reform[data$t2==1] + ps_reform[data$t2==1] + qr_reform[data$t2==1] + tardecreases[data$t2==1] + com.t2 + factor(region[data$t2==1]))
m1_reg_2_t2 <- lm(reg_likeli2[data$t2==1] ~ tarmeasure[data$t2==1] + fdi_reform[data$t2==1] + licensing_reform[data$t2==1] + ps_reform[data$t2==1] + qr_reform[data$t2==1] + tardecreases[data$t2==1] + reg.t2 + factor(region[data$t2==1]))

fv_bjp_t2_ldv <- m1_bjp_2_t2$fitted.values
fv_com_t2_ldv <- m1_com_2_t2$fitted.values
fv_reg_t2_ldv <- m1_reg_2_t2$fitted.values
fv_inc_t2_ldv <- 0 - (m1_bjp_2_t2$fitted.values + m1_com_2_t2$fitted.values + m1_reg_2_t2$fitted.values)

untran_fv_bjp_t2_ldv <- exp(fv_bjp_t2_ldv)/(1+exp(fv_bjp_t2_ldv)+exp(fv_com_t2_ldv)+exp(fv_reg_t2_ldv))
untran_fv_com_t2_ldv <- exp(fv_com_t2_ldv)/(1+exp(fv_bjp_t2_ldv)+exp(fv_com_t2_ldv)+exp(fv_reg_t2_ldv))
untran_fv_reg_t2_ldv <- exp(fv_reg_t2_ldv)/(1+exp(fv_bjp_t2_ldv)+exp(fv_com_t2_ldv)+exp(fv_reg_t2_ldv))
untran_fv_inc_t2_ldv <- 1-(untran_fv_bjp_t2_ldv + untran_fv_com_t2_ldv + untran_fv_reg_t2_ldv)

exp_bjp_t2_ldv <- (untran_fv_bjp_t2_ldv*4)-1
exp_com_t2_ldv <- (untran_fv_com_t2_ldv*4)-1
exp_reg_t2_ldv <- (untran_fv_reg_t2_ldv*4)-1
exp_inc_t2_ldv <- (untran_fv_inc_t2_ldv*4)-1

exp_bjpreg_t2_ldv <- exp_bjp_t2_ldv+exp_reg_t2_ldv

hex_fvs_t2_ldv <- terncoord(cbind(exp_bjpreg_t2_ldv,exp_inc_t2_ldv,exp_com_t2_ldv))

##PLOT 2
pdf(file= "Hextime_ldv.pdf", width = 5, height = 5, family = "Helvetica", pointsize = 10)
par(mfrow=c(2,2),mar=c(1,1,1,1))
plot(hex_fvs_t2_ldv[,1],hex_fvs_t2_ldv[,2],xlim=c(-1.2,1.2),ylim=c(-1.2,1.2),pch=3,col="blue",cex=0.5,xlab="",ylab="",axes=FALSE)
points(x=l1[,1],y=l1[,2])
points(x=l2[,1],y=l2[,2])
points(x=l3[,1],y=l3[,2])
points(x=l4[,1],y=l4[,2])
points(x=l5[,1],y=l5[,2])
points(x=l6[,1],y=l6[,2])
polygon(x=c(-1.29,1.29,1.29,-1.29),y=c(1.29,1.29,-1.29,-1.29))

#Full hexagon contours - 20%iles
slope <- sqrt(3)
twentyunit <- sqrt(3)/5
for (i in -5:5){
  abline(i*twentyunit,slope,lty=2, col="darkgray")
  abline(i*twentyunit,-slope,lty=2, col="darkgray")
  abline(i*twentyunit/2,0,lty=2, col="darkgray")
}
abline(0*twentyunit,slope,lty=1,lwd=2)
abline(0*twentyunit,-slope,lty=1,lwd=2)
abline(0*twentyunit/2,0,lty=1,lwd=2)
lines(x=c(l1[,1],l2[,1]),y=c(l3[,2],l3[,2]),lty=1,lwd=2,col="green")
lines(x=c(l1[,1],l2[,1]),y=c(l4[,2],l4[,2]),lty=1,lwd=2,col="red")
lines(x=c(l1[,1],l5[,1]),y=c(l4[,2],l5[,2]),lty=1,lwd=2,col="green")
lines(x=c(l2[,1],l6[,1]),y=c(l1[,2],l5[,2]),lty=1,lwd=2,col="red")
lines(x=c(l6[,1],l3[,1]),y=c(l5[,2],l4[,2]),lty=1,lwd=2,col="green")
lines(x=c(l4[,1],l5[,1]),y=c(l3[,2],l5[,2]),lty=1,lwd=2,col="red")

#Interior hexagon contours
#Labels
text(0,1.1,"+ BJP/Regional")
text(-0.85,-0.5,"+ INC",srt=-57.5)
text(0.85,-0.5,"+ Far Left",srt=57.5)
text(0,-1.1,"- BJP/Regional")
text(+0.85,+0.5,"- INC",srt=-57.5)
text(-0.90,+0.5,"- Far Left",srt=57.5)
text(-1,+1.1,"1996-1998")

#####By time period #t=3#Time FEs removed from models.
m1_bjp_2_t3 <- lm(bjp_likeli2[data$t3==1] ~ tarmeasure[data$t3==1] + fdi_reform[data$t3==1] + licensing_reform[data$t3==1] + ps_reform[data$t3==1] + qr_reform[data$t3==1] + tardecreases[data$t3==1] + bjp.t3 + factor(region[data$t3==1]))
m1_com_2_t3 <- lm(com_likeli2[data$t3==1] ~ tarmeasure[data$t3==1] + fdi_reform[data$t3==1] + licensing_reform[data$t3==1] + ps_reform[data$t3==1] + qr_reform[data$t3==1] + tardecreases[data$t3==1] + com.t3 + factor(region[data$t3==1]))
m1_reg_2_t3 <- lm(reg_likeli2[data$t3==1] ~ tarmeasure[data$t3==1] + fdi_reform[data$t3==1] + licensing_reform[data$t3==1] + ps_reform[data$t3==1] + qr_reform[data$t3==1] + tardecreases[data$t3==1] + reg.t3 + factor(region[data$t3==1]))

fv_bjp_t3_ldv <- m1_bjp_2_t3$fitted.values
fv_com_t3_ldv <- m1_com_2_t3$fitted.values
fv_reg_t3_ldv <- m1_reg_2_t3$fitted.values
fv_inc_t3_ldv <- 0 - (m1_bjp_2_t3$fitted.values + m1_com_2_t3$fitted.values + m1_reg_2_t3$fitted.values)

untran_fv_bjp_t3_ldv <- exp(fv_bjp_t3_ldv)/(1+exp(fv_bjp_t3_ldv)+exp(fv_com_t3_ldv)+exp(fv_reg_t3_ldv))
untran_fv_com_t3_ldv <- exp(fv_com_t3_ldv)/(1+exp(fv_bjp_t3_ldv)+exp(fv_com_t3_ldv)+exp(fv_reg_t3_ldv))
untran_fv_reg_t3_ldv <- exp(fv_reg_t3_ldv)/(1+exp(fv_bjp_t3_ldv)+exp(fv_com_t3_ldv)+exp(fv_reg_t3_ldv))
untran_fv_inc_t3_ldv <- 1-(untran_fv_bjp_t3_ldv + untran_fv_com_t3_ldv + untran_fv_reg_t3_ldv)

exp_bjp_t3_ldv <- (untran_fv_bjp_t3_ldv*4)-1
exp_com_t3_ldv <- (untran_fv_com_t3_ldv*4)-1
exp_reg_t3_ldv <- (untran_fv_reg_t3_ldv*4)-1
exp_inc_t3_ldv <- (untran_fv_inc_t3_ldv*4)-1

exp_bjpreg_t3_ldv <- exp_bjp_t3_ldv+exp_reg_t3_ldv

hex_fvs_t3_ldv <- terncoord(cbind(exp_bjpreg_t3_ldv,exp_inc_t3_ldv,exp_com_t3_ldv))

plot(hex_fvs_t3_ldv[,1],hex_fvs_t3_ldv[,2],xlim=c(-1.2,1.2),ylim=c(-1.2,1.2),pch=3,col="blue",cex=0.5,xlab="",ylab="",axes=FALSE)
points(x=l1[,1],y=l1[,2])
points(x=l2[,1],y=l2[,2])
points(x=l3[,1],y=l3[,2])
points(x=l4[,1],y=l4[,2])
points(x=l5[,1],y=l5[,2])
points(x=l6[,1],y=l6[,2])
polygon(x=c(-1.29,1.29,1.29,-1.29),y=c(1.29,1.29,-1.29,-1.29))

#Full hexagon contours - 20%iles
slope <- sqrt(3)
twentyunit <- sqrt(3)/5
for (i in -5:5){
  abline(i*twentyunit,slope,lty=2, col="darkgray")
  abline(i*twentyunit,-slope,lty=2, col="darkgray")
  abline(i*twentyunit/2,0,lty=2, col="darkgray")
}
abline(0*twentyunit,slope,lty=1,lwd=2)
abline(0*twentyunit,-slope,lty=1,lwd=2)
abline(0*twentyunit/2,0,lty=1,lwd=2)
lines(x=c(l1[,1],l2[,1]),y=c(l3[,2],l3[,2]),lty=1,lwd=2,col="green")
lines(x=c(l1[,1],l2[,1]),y=c(l4[,2],l4[,2]),lty=1,lwd=2,col="red")
lines(x=c(l1[,1],l5[,1]),y=c(l4[,2],l5[,2]),lty=1,lwd=2,col="green")
lines(x=c(l2[,1],l6[,1]),y=c(l1[,2],l5[,2]),lty=1,lwd=2,col="red")
lines(x=c(l6[,1],l3[,1]),y=c(l5[,2],l4[,2]),lty=1,lwd=2,col="green")
lines(x=c(l4[,1],l5[,1]),y=c(l3[,2],l5[,2]),lty=1,lwd=2,col="red")

#Interior hexagon contours
#Labels
text(0,1.1,"+ BJP/Regional")
text(-0.85,-0.5,"+ INC",srt=-57.5)
text(0.85,-0.5,"+ Far Left",srt=57.5)
text(0,-1.1,"- BJP/Regional")
text(+0.85,+0.5,"- INC",srt=-57.5)
text(-0.90,+0.5,"- Far Left",srt=57.5)
text(-1,+1.1,"1998-1999")

#####By time period #t=4#Time FEs removed from models.
m1_bjp_2_t4 <- lm(bjp_likeli2[data$t4==1] ~ tarmeasure[data$t4==1] + fdi_reform[data$t4==1] + licensing_reform[data$t4==1] + ps_reform[data$t4==1] + qr_reform[data$t4==1] + tardecreases[data$t4==1] + bjp.t4 + factor(region[data$t4==1]))
m1_com_2_t4 <- lm(com_likeli2[data$t4==1] ~ tarmeasure[data$t4==1] + fdi_reform[data$t4==1] + licensing_reform[data$t4==1] + ps_reform[data$t4==1] + qr_reform[data$t4==1] + tardecreases[data$t4==1] + com.t4 + factor(region[data$t4==1]))
m1_reg_2_t4 <- lm(reg_likeli2[data$t4==1] ~ tarmeasure[data$t4==1] + fdi_reform[data$t4==1] + licensing_reform[data$t4==1] + ps_reform[data$t4==1] + qr_reform[data$t4==1] + tardecreases[data$t4==1] + reg.t4 + factor(region[data$t4==1]))

fv_bjp_t4_ldv <- m1_bjp_2_t4$fitted.values
fv_com_t4_ldv <- m1_com_2_t4$fitted.values
fv_reg_t4_ldv <- m1_reg_2_t4$fitted.values
fv_inc_t4_ldv <- 0 - (m1_bjp_2_t4$fitted.values + m1_com_2_t4$fitted.values + m1_reg_2_t4$fitted.values)

untran_fv_bjp_t4_ldv <- exp(fv_bjp_t4_ldv)/(1+exp(fv_bjp_t4_ldv)+exp(fv_com_t4_ldv)+exp(fv_reg_t4_ldv))
untran_fv_com_t4_ldv <- exp(fv_com_t4_ldv)/(1+exp(fv_bjp_t4_ldv)+exp(fv_com_t4_ldv)+exp(fv_reg_t4_ldv))
untran_fv_reg_t4_ldv <- exp(fv_reg_t4_ldv)/(1+exp(fv_bjp_t4_ldv)+exp(fv_com_t4_ldv)+exp(fv_reg_t4_ldv))
untran_fv_inc_t4_ldv <- 1-(untran_fv_bjp_t4_ldv + untran_fv_com_t4_ldv + untran_fv_reg_t4_ldv)

exp_bjp_t4_ldv <- (untran_fv_bjp_t4_ldv*4)-1
exp_com_t4_ldv <- (untran_fv_com_t4_ldv*4)-1
exp_reg_t4_ldv <- (untran_fv_reg_t4_ldv*4)-1
exp_inc_t4_ldv <- (untran_fv_inc_t4_ldv*4)-1

exp_bjpreg_t4_ldv <- exp_bjp_t4_ldv+exp_reg_t4_ldv

hex_fvs_t4_ldv <- terncoord(cbind(exp_bjpreg_t4_ldv,exp_inc_t4_ldv,exp_com_t4_ldv))

plot(hex_fvs_t4_ldv[,1],hex_fvs_t4_ldv[,2],xlim=c(-1.2,1.2),ylim=c(-1.2,1.2),pch=3,col="blue",cex=0.5,xlab="",ylab="",axes=FALSE)
points(x=l1[,1],y=l1[,2])
points(x=l2[,1],y=l2[,2])
points(x=l3[,1],y=l3[,2])
points(x=l4[,1],y=l4[,2])
points(x=l5[,1],y=l5[,2])
points(x=l6[,1],y=l6[,2])
polygon(x=c(-1.29,1.29,1.29,-1.29),y=c(1.29,1.29,-1.29,-1.29))

#Full hexagon contours - 20%iles
slope <- sqrt(3)
twentyunit <- sqrt(3)/5
for (i in -5:5){
  abline(i*twentyunit,slope,lty=2, col="darkgray")
  abline(i*twentyunit,-slope,lty=2, col="darkgray")
  abline(i*twentyunit/2,0,lty=2, col="darkgray")
}
abline(0*twentyunit,slope,lty=1,lwd=2)
abline(0*twentyunit,-slope,lty=1,lwd=2)
abline(0*twentyunit/2,0,lty=1,lwd=2)
lines(x=c(l1[,1],l2[,1]),y=c(l3[,2],l3[,2]),lty=1,lwd=2,col="green")
lines(x=c(l1[,1],l2[,1]),y=c(l4[,2],l4[,2]),lty=1,lwd=2,col="red")
lines(x=c(l1[,1],l5[,1]),y=c(l4[,2],l5[,2]),lty=1,lwd=2,col="green")
lines(x=c(l2[,1],l6[,1]),y=c(l1[,2],l5[,2]),lty=1,lwd=2,col="red")
lines(x=c(l6[,1],l3[,1]),y=c(l5[,2],l4[,2]),lty=1,lwd=2,col="green")
lines(x=c(l4[,1],l5[,1]),y=c(l3[,2],l5[,2]),lty=1,lwd=2,col="red")

#Interior hexagon contours
#Labels
text(0,1.1,"+ BJP/Regional")
text(-0.85,-0.5,"+ INC",srt=-57.5)
text(0.85,-0.5,"+ Far Left",srt=57.5)
text(0,-1.1,"- BJP/Regional")
text(+0.85,+0.5,"- INC",srt=-57.5)
text(-0.90,+0.5,"- Far Left",srt=57.5)
text(-1,+1.1,"1999-2004")
dev.off()



#############################################################

##FIGURE 8: First Differences
bjpreg_fd_ldv <- fd2_bjp_2 + fd2_reg_2
inc_fd_ldv <- fd2_inc_2
com_fd_ldv <- fd2_com_2
fd_ldv <- cbind(bjpreg_fd_ldv,inc_fd_ldv,com_fd_ldv)

hex_fd_ldv <- terncoord(fd_ldv)

par(mfrow=c(1,1))
pdf(file= "Hexfd_ldv.pdf", width = 5, height = 5, family = "Helvetica", pointsize = 10)
plot(hex_fd_ldv[,1],hex_fd_ldv[,2],xlim=c(-1.2,1.2),ylim=c(-1.2,1.2),pch=3,col="blue",cex=0.5,xlab="",ylab="",axes=FALSE)
points(x=l1[,1],y=l1[,2])
points(x=l2[,1],y=l2[,2])
points(x=l3[,1],y=l3[,2])
points(x=l4[,1],y=l4[,2])
points(x=l5[,1],y=l5[,2])
points(x=l6[,1],y=l6[,2])
polygon(x=c(-1.29,1.29,1.29,-1.29),y=c(1.29,1.29,-1.29,-1.29))


#Full hexagon contours - 20%iles
slope <- sqrt(3)
twentyunit <- sqrt(3)/5
for (i in -5:5){
  abline(i*twentyunit,slope,lty=2, col="darkgray")
  abline(i*twentyunit,-slope,lty=2, col="darkgray")
  abline(i*twentyunit/2,0,lty=2, col="darkgray")
}
abline(0*twentyunit,slope,lty=1,lwd=2)
abline(0*twentyunit,-slope,lty=1,lwd=2)
abline(0*twentyunit/2,0,lty=1,lwd=2)
lines(x=c(l1[,1],l2[,1]),y=c(l3[,2],l3[,2]),lty=1,lwd=2,col="green")
lines(x=c(l1[,1],l2[,1]),y=c(l4[,2],l4[,2]),lty=1,lwd=2,col="red")
lines(x=c(l1[,1],l5[,1]),y=c(l4[,2],l5[,2]),lty=1,lwd=2,col="green")
lines(x=c(l2[,1],l6[,1]),y=c(l1[,2],l5[,2]),lty=1,lwd=2,col="red")
lines(x=c(l6[,1],l3[,1]),y=c(l5[,2],l4[,2]),lty=1,lwd=2,col="green")
lines(x=c(l4[,1],l5[,1]),y=c(l3[,2],l5[,2]),lty=1,lwd=2,col="red")

#Interior hexagon contours
#TBC

#Labels
text(0,1.1,"+ BJP/Regional")
text(-0.85,-0.5,"+ INC",srt=-57.5)
text(0.85,-0.5,"+ Far Left",srt=57.5)
text(0,-1.1,"- BJP/Regional")
text(+0.85,+0.5,"- INC",srt=-57.5)
text(-0.90,+0.5,"- Far Left",srt=57.5)
dev.off()

#############################################################

## FIGURE 9: First Differences: Group Far Let
bjp_fd_ldv2<-fd2_bjp_2
reg_fd_ldv2<-fd2_reg_2
left_fd_ldv2<-fd2_inc_2 + fd2_com_2
fd_ldv2<-cbind(bjp_fd_ldv2, left_fd_ldv2, reg_fd_ldv2)

hex_fd_ldv2 <- terncoord(fd_ldv2)

par(mfrow=c(1,1))
pdf(file= "Hexfd_left.pdf", width = 5, height = 5, family = "Helvetica", pointsize = 10)
plot(hex_fd_ldv2[,1],hex_fd_ldv2[,2],xlim=c(-1.2,1.2),ylim=c(-1.2,1.2),pch=3,col="blue",cex=0.5,xlab="",ylab="",axes=FALSE)
points(x=l1[,1],y=l1[,2])
points(x=l2[,1],y=l2[,2])
points(x=l3[,1],y=l3[,2])
points(x=l4[,1],y=l4[,2])
points(x=l5[,1],y=l5[,2])
points(x=l6[,1],y=l6[,2])
polygon(x=c(-1.29,1.29,1.29,-1.29),y=c(1.29,1.29,-1.29,-1.29))


#Full hexagon contours - 20%iles
slope <- sqrt(3)
twentyunit <- sqrt(3)/5
for (i in -5:5){
  abline(i*twentyunit,slope,lty=2, col="darkgray")
  abline(i*twentyunit,-slope,lty=2, col="darkgray")
  abline(i*twentyunit/2,0,lty=2, col="darkgray")
}
abline(0*twentyunit,slope,lty=1,lwd=2)
abline(0*twentyunit,-slope,lty=1,lwd=2)
abline(0*twentyunit/2,0,lty=1,lwd=2)
lines(x=c(l1[,1],l2[,1]),y=c(l3[,2],l3[,2]),lty=1,lwd=2,col="green")
lines(x=c(l1[,1],l2[,1]),y=c(l4[,2],l4[,2]),lty=1,lwd=2,col="red")
lines(x=c(l1[,1],l5[,1]),y=c(l4[,2],l5[,2]),lty=1,lwd=2,col="green")
lines(x=c(l2[,1],l6[,1]),y=c(l1[,2],l5[,2]),lty=1,lwd=2,col="red")
lines(x=c(l6[,1],l3[,1]),y=c(l5[,2],l4[,2]),lty=1,lwd=2,col="green")
lines(x=c(l4[,1],l5[,1]),y=c(l3[,2],l5[,2]),lty=1,lwd=2,col="red")

#Labels
text(0,1.1,"+ BJP")
text(-0.85,-0.5,"+ Left",srt=-57.5)
text(0.85,-0.5,"+ Regional",srt=57.5)
text(0,-1.1,"- BJP")
text(+0.85,+0.5,"- Left",srt=-57.5)
text(-0.90,+0.5,"- Regional",srt=57.5)
dev.off()






#############################################################
##---------------------Appendix A--------------------------##
##---------Replicating Tandon (2012) Tables 5-7------------##
#############################################################

##TABLE 3: Replicating Table 5- Support for the far left and congress in response to the tariff reform

##Model 1
##Since the data for this model is panel in structure but does not include group fixed effects, model="pooling" is specified to just run a panel regression model
lm5.1<-plm(left~tarmeasure + factor(time), model="pooling", data=data.plm)
lm5.1
##Standard errors are calculated as normal
se5.1<-vcovHC(lm5.1, method="arellano", type="HC1", cluster="group")*(m/(m-1))
se5.1<-sqrt(diag(se5.1))
se5.1

##Model 2
##This model returns to the method used in tables 3 and 4, with method="arellano"
lm5.2<-plm(left~tarmeasure + tardecreases + fdi_reform + ps_reform + qr_reform + licensing_reform + factor(time), model="within", data=data.plm)
lm5.2
se5.2<-vcovHC(lm5.2, method="arellano", type="HC1", cluster="group")*(m/(m-1))
se5.2<-sqrt(diag(se5.2))
se5.2

##Model 3
iv5.3<-ivreg(left~tarmeasure + tardecreases + fdi_reform + ps_reform + qr_reform + licensing_reform + factor(time) + factor(region) | tarmeasure_instrument + tardecreases + fdi_reform + ps_reform + qr_reform + licensing_reform + factor(time) + factor(region), data=data.plm, model=T)
summary(iv5.3)

##Model 4
##Using the tobit function in the AER package, this model calculates the tobit with limits of -1 and 1
tobit5.4<- tobit(left~tarmeasure + tardecreases + fdi_reform + ps_reform + qr_reform + licensing_reform + factor(time), left=-1, right=1, data=data.plm)
summary(tobit5.4)

#############################################################

##TABLE 4: Replicating Table 6- Cross-sectional voting response to reforms between 1991 and 1996
##Cross Sectional analyses using Only Time 1

##Subset Data for each time frame
time1<-data[t2==0 & t3==0 & t4==0,] 

##Model 1
##Since this model is cross-sectional, we return to the traditional lm function and the corresponding clustered standard error function
lm6.1<-lm(left~tarmeasure, data=time1)
lm6.1
clx(lm6.1, 1, time1$region)

##Model 2
lm6.2<-lm(left~tarmeasure + tardecreases + fdi_reform + ps_reform + qr_reform + licensing_reform + coast + rural + married + log_landown + illiterate + poor, data=time1)
lm6.2
clx(lm6.2, 1, time1$region)

##Model 3
iv6.3<-ivreg(left~tarmeasure + tardecreases + fdi_reform + ps_reform + qr_reform + licensing_reform + coast + rural + married + log_landown + illiterate + poor | tarmeasure_instrument + tardecreases + fdi_reform + ps_reform + qr_reform + licensing_reform + coast + rural + married + log_landown + illiterate + poor, model=T, data=time1)
summary(iv6.3)
clx(iv6.3, 1, time1$region)

##Model 4
lm6.4<-lm(bjp~tarmeasure, data=time1)
lm6.4
clx(lm6.4, 1, time1$region)

##Model 5
lm6.5<-lm(bjp~tarmeasure + tardecreases + fdi_reform + ps_reform + qr_reform + licensing_reform + coast + rural + married + log_landown + illiterate + poor, data=time1)
lm6.5
clx(lm6.5, 1, time1$region)

##Model 6
iv6.6<-ivreg(bjp~tarmeasure + tardecreases + fdi_reform + ps_reform + qr_reform + licensing_reform + coast + rural + married + log_landown + illiterate + poor | tarmeasure_instrument + tardecreases + fdi_reform + ps_reform + qr_reform + licensing_reform + coast + rural + married + log_landown + illiterate + poor, model=T, data=time1)
summary(lm6.6)
clx(iv6.6, 1, time1$region)

#############################################################

##TABLE 5: Replicating Table 7- Vote changes prior to the reform

##Model 1
##These models return to the orginal baseline data and therefore employ a panel structure
lm7.1<-plm(leftvote_2~tarmeasure + tardecreases + fdi_reform + ps_reform + qr_reform + licensing_reform + factor(time), model="within", data=data.plm)
lm7.1
se7.1<-vcovHC(lm7.1, method="arellano", type="HC1", cluster="group")*(m/(m-1))
se7.1<-sqrt(diag(se7.1))
se7.1

##Model 2
iv7.2<-ivreg(leftvote_2~tarmeasure + tardecreases + fdi_reform + ps_reform + qr_reform + licensing_reform + factor(time) + factor(region) | tarmeasure_instrument + tardecreases + fdi_reform + ps_reform + qr_reform + licensing_reform + factor(time) + factor(region), data=data.plm)
summary(iv7.2)

##Model 3
lm7.3<-plm(leftvote_1~tarmeasure + tardecreases + fdi_reform + ps_reform + qr_reform + licensing_reform + factor(time), model="within", data=data.plm)
lm7.3
se7.3<-vcovHC(lm7.3, method="arellano", type="HC1", cluster="group")*(m/(m-1))
se7.3<-sqrt(diag(se7.3))
se7.3

##Model 4
iv7.4<-ivreg(leftvote_1~tarmeasure + tardecreases + fdi_reform + ps_reform + qr_reform + licensing_reform + factor(time) + factor(region) | tarmeasure_instrument + tardecreases + fdi_reform + ps_reform + qr_reform + licensing_reform + factor(time) + factor(region), data=data.plm)
summary(iv7.4)

##Model 5
lm7.5<-plm(left~tarmeasure + tardecreases_future + tardecreases + fdi_reform + ps_reform + qr_reform + licensing_reform + factor(time), model="within", data=data.plm)
lm7.5
se7.5<-vcovHC(lm7.5, method="arellano", type="HC1", cluster="group")*(m/(m-1))
se7.5<-sqrt(diag(se7.5))
se7.5

##Model 6
iv7.6<-ivreg(left~tarmeasure + tardecreases_future +  tardecreases + fdi_reform + ps_reform + qr_reform + licensing_reform + factor(time) + factor(region) | tarmeasure_instrument + tardecreases_future + tardecreases + fdi_reform + ps_reform + qr_reform + licensing_reform + factor(time) + factor(region), data=data.plm)
summary(iv7.6)


#############################################################
##---------------------Appendix D--------------------------##
##------Separating the INC and Far Left: Tables 5-7--------##
#############################################################


##TABLE 6: INC v. FAR LEFT Table 5- Support for the far left and congress in response to the tariff reform

##INC
##Model 1
##Since the data for this model is panel in structure but does not include group fixed effects, model="pooling" is specified to just run a panel regression model
lm5.1a<-plm(inc~tarmeasure + factor(time), model="pooling", data=data.plm)
lm5.1a
##Standard errors are calculated as normal
se5.1a<-vcovHC(lm5.1a, method="arellano", type="HC1", cluster="group")*(m/(m-1))
se5.1a<-sqrt(diag(se5.1a))
se5.1a

##Model 2
##This model returns to the method used in tables 3 and 4, with method="arellano"
lm5.2a<-plm(inc~tarmeasure + tardecreases + fdi_reform + ps_reform + qr_reform + licensing_reform + factor(time), model="within", data=data.plm)
lm5.2a
se5.2a<-vcovHC(lm5.2a, method="arellano", type="HC1", cluster="group")*(m/(m-1))
se5.2a<-sqrt(diag(se5.2a))
se5.2a

##Model 3
iv5.3a<-ivreg(inc~tarmeasure + tardecreases + fdi_reform + ps_reform + qr_reform + licensing_reform + factor(time) + factor(region) | tarmeasure_instrument + tardecreases + fdi_reform + ps_reform + qr_reform + licensing_reform + factor(time) + factor(region), data=data.plm, model=T)
summary(iv5.3a)

##Model 4
##Using the tobit function in the AER package, this model calculates the tobit with limits of -1 and 1
tobit5.4a<- tobit(inc~tarmeasure + tardecreases + fdi_reform + ps_reform + qr_reform + licensing_reform + factor(time), left=-1, right=1, data=data.plm)
summary(tobit5.4a)



##FAR LEFT
##Model 1
##Since the data for this model is panel in structure but does not include group fixed effects, model="pooling" is specified to just run a panel regression model
lm5.1b<-plm(com~tarmeasure + factor(time), model="pooling", data=data.plm)
lm5.1b
##Standard errors are calculated as normal
se5.1b<-vcovHC(lm5.1b, method="arellano", type="HC1", cluster="group")*(m/(m-1))
se5.1b<-sqrt(diag(se5.1b))
se5.1b

##Model 2
##This model returns to the method used in tables 3 and 4, with method="arellano"
lm5.2b<-plm(com~tarmeasure + tardecreases + fdi_reform + ps_reform + qr_reform + licensing_reform + factor(time), model="within", data=data.plm)
lm5.2b
se5.2b<-vcovHC(lm5.2b, method="arellano", type="HC1", cluster="group")*(m/(m-1))
se5.2b<-sqrt(diag(se5.2b))
se5.2b

##Model 3
iv5.3b<-ivreg(com~tarmeasure + tardecreases + fdi_reform + ps_reform + qr_reform + licensing_reform + factor(time) + factor(region) | tarmeasure_instrument + tardecreases + fdi_reform + ps_reform + qr_reform + licensing_reform + factor(time) + factor(region), data=data.plm, model=T)
summary(iv5.3b)

##Model 4
##Using the tobit function in the AER package, this model calculates the tobit with limits of -1 and 1
tobit5.4b<- tobit(com~tarmeasure + tardecreases + fdi_reform + ps_reform + qr_reform + licensing_reform + factor(time), left=-1, right=1, data=data.plm)
summary(tobit5.4b)

#############################################################

##TABLE 7: INC v. FAR LEFT Table 6- Cross-sectional voting response to reforms between 1991 and 1996
##Cross Sectional analyses using Only Time 1

##INC
##Model 1
##Since this model is cross-sectional, we return to the traditional lm function and the corresponding clustered standard error function
lm6.1a<-lm(inc~tarmeasure, data=time1)
lm6.1a
clx(lm6.1a, 1, time1$region)

##Model 2
lm6.2a<-lm(inc~tarmeasure + tardecreases + fdi_reform + ps_reform + qr_reform + licensing_reform + coast + rural + married + log_landown + illiterate + poor, data=time1)
lm6.2a
clx(lm6.2a, 1, time1$region)

##Model 3
iv6.3a<-ivreg(inc~tarmeasure + tardecreases + fdi_reform + ps_reform + qr_reform + licensing_reform + coast + rural + married + log_landown + illiterate + poor | tarmeasure_instrument + tardecreases + fdi_reform + ps_reform + qr_reform + licensing_reform + coast + rural + married + log_landown + illiterate + poor, model=T, data=time1)
summary(iv6.3a)
clx(iv6.3a, 1, time1$region)



##FAR LEFT
##Model 1
##Since this model is cross-sectional, we return to the traditional lm function and the corresponding clustered standard error function
lm6.1b-lm(com~tarmeasure, data=time1)
lm6.1b
clx(lm6.1b, 1, time1$region)

##Model 2
lm6.2b<-lm(com~tarmeasure + tardecreases + fdi_reform + ps_reform + qr_reform + licensing_reform + coast + rural + married + log_landown + illiterate + poor, data=time1)
lm6.2b
clx(lm6.2b, 1, time1$region)

##Model 3
iv6.3b<-ivreg(com~tarmeasure + tardecreases + fdi_reform + ps_reform + qr_reform + licensing_reform + coast + rural + married + log_landown + illiterate + poor | tarmeasure_instrument + tardecreases + fdi_reform + ps_reform + qr_reform + licensing_reform + coast + rural + married + log_landown + illiterate + poor, model=T, data=time1)
summary(iv6.3b)
clx(iv6.3b, 1, time1$region)


#############################################################

##TABLE 8: INC v. FAR LEFT Table 7- Vote changes prior to the reform

##Generate Lags
data.plm$inc_2<-lag(data.plm$inc,2)
data.plm$com_2<-lag(data.plm$com,2)


##INC
##Model 1
##These models return to the orginal baseline data and therefore employ a panel structure
lm7.1a<-plm(inc_2~tarmeasure + tardecreases + fdi_reform + ps_reform + qr_reform + licensing_reform + factor(time), model="within", data=data.plm)
lm7.1a
se7.1a<-vcovHC(lm7.1a, method="arellano", type="HC1", cluster="group")*(m/(m-1))
se7.1a<-sqrt(diag(se7.1a))
se7.1a

##Model 3
lm7.3a<-plm(inc_1~tarmeasure + tardecreases + fdi_reform + ps_reform + qr_reform + licensing_reform + factor(time), model="within", data=data.plm)
lm7.3a
se7.3a<-vcovHC(lm7.3a, method="arellano", type="HC1", cluster="group")*(m/(m-1))
se7.3a<-sqrt(diag(se7.3a))
se7.3a

##Model 5
lm7.5a<-plm(inc~tarmeasure + tardecreases_future + tardecreases + fdi_reform + ps_reform + qr_reform + licensing_reform + factor(time), model="within", data=data.plm)
lm7.5a
se7.5a<-vcovHC(lm7.5a, method="arellano", type="HC1", cluster="group")*(m/(m-1))
se7.5a<-sqrt(diag(se7.5a))
se7.5a


##FAR LEFT
##Model 1
##These models return to the orginal baseline data and therefore employ a panel structure
lm7.1b<-plm(com_2~tarmeasure + tardecreases + fdi_reform + ps_reform + qr_reform + licensing_reform + factor(time), model="within", data=data.plm)
lm7.1b
se7.1b<-vcovHC(lm7.1b, method="arellano", type="HC1", cluster="group")*(m/(m-1))
se7.1b<-sqrt(diag(se7.1b))
se7.1b

##Model 3
lm7.3b<-plm(com_1~tarmeasure + tardecreases + fdi_reform + ps_reform + qr_reform + licensing_reform + factor(time), model="within", data=data.plm)
lm7.3b
se7.3b<-vcovHC(lm7.3b, method="arellano", type="HC1", cluster="group")*(m/(m-1))
se7.3b<-sqrt(diag(se7.3b))
se7.3b

##Model 5
lm7.5b<-plm(com~tarmeasure + tardecreases_future + tardecreases + fdi_reform + ps_reform + qr_reform + licensing_reform + factor(time), model="within", data=data.plm)
lm7.5b
se7.5b<-vcovHC(lm7.5b, method="arellano", type="HC1", cluster="group")*(m/(m-1))
se7.5b<-sqrt(diag(se7.5b))
se7.5b




